CONTROL OF NONLINEAR SYSTEMS 
REPRESENTED IN QUASILINEAR FORM 
Final Report 

NASA Grant NAG-1-126 
Josef A. Coetsee 


December 1993 




CONTROL OF NONLINEAR SYSTEMS 
REPRESENTED IN QUASILINEAR 

FORM 

by 

Josef Adriaan Coetsee 

B. Sc. (Hons) Elek. Rand Afrikaans University, 1979 
S.M. Aeronautics and Astronautics Massachusetts Institute of Technology. 1985 

SUBMITTED TO THE DEPARTMENT OF 
AERONAUTICS AND ASTRONAUTICS IN PARTIAL FULFILLMENT OF THE 
REQUIREMENTS FOR THE DEGREE OF 

Doctor of Philosophy 
in Aeronautics and Astronautics 

at the 

Massachusetts Institute of Technology 

February 1994 

© Massachusetts Institute of Technology 1994. All rights reserved. 


Signature of Author. 


Certified bv 


Certified bv 


Certified bv 


Certified bv 






Department ot Aeronautics And Astronautics 

tti 

Prof. W.E. Vander Velde, Thesis Supervisor 
Department of Aeronautics and Astronautics 

//: / 


s/ 

Dr. H. Alexander, Thesis Committee Member 
Senioi^Engineer, Orbital Sciery^ Cm^poration 


. 

Prof. J-J.E. Slotine, Thesis Committee Member 
Department of Mechanical Engineering 

1 }/ 

Dr. A. Von Flotow, Thesis Committee Member 
President, Hood Technology Corporation 


Accepted by 


Prof. Harold Y. Wachman. Chairman 
Department Graduate Committee 



CONTROL OF NONLINEAR SYSTEMS REPRESENTED 

IN QUASILINEAR FORM 

by 

Josef Adriaan Coetsee 

Submitted to the Department of Aeronautics And Astronautics 
on December 17, 1993, in partial fulfillment of the 
requirements for the degree of 
Doctor of Philosophy in Aeronautics and Astronautics 


Abstract 


Methods to synthesize controllers for nonlinear systems are developed by exploiting 
the fact that under mild differentiability conditions, systems of the form: 

x = f(x) + G(x.)u 

can be represented in quasilinear form, viz: 

x = A(x)x 4- 5(x) u 

Two classes of control methods are investigated: 

• zero-look-ahead control, where the control input depends only on the current val- 
ues of A(x), 5(x). For this case the control input is computed by continuously 
solving a matrix Ricatti equation as the system progresses along a trajectory. 

• controllers with look-ahead, where the control input depends on the future be- 
havior of A(x),J3(x). These controllers use the similarity between quasilinear 
systems, and linear time varying systems to find approximate solutions to op- 
timal control type problems. 

The methods that are developed are not guaranteed to be globally stable. However in 
simulation studies they were found to be useful alternatives for synthesizing control 
laws for a general class of nonlinear systems. 


Prof. W.E. Vander Velde, Thesis Supervisor 
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Chapter 1 


Introduction 


Control engineering has played an important role in the development of industrial 
society. One of the earliest applications of control devices can be found in Watt's 
use of flyball governors to control the speed of steam engines [42]. Subsequently 
the discipline has undergone a great deal of theoretical development and has found 
practical application in a wide variety of systems such as the control of chemical 
processes, robotic manipulators, flight vehicles etc. 

To develop control systems the following paradigm is often used: 

1. It is assumed that there is a requirement for a system, also referred to as the 
plant, to behave in accordance with certain specifications. For example the 
requirement may be for a motor (the plant) to operate at a certain speed. 

2. When the system s natural behavior does not in itself meet the requirements, it 
may be necessary to use some control inputs to induce the desired behavior. As 
an example some “regulator” may be required to ensure that a motor operates 
at the required speed. 

3. The next step is to obtain an abstract dynamical model for the system. Usually 
it is assumed that the system's behavior is adequately described by a set of 
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ordinary differential equations. The differential equations are often derived on 
the basis of some physical principles, but other approaches are possible. For 
example it can simply be assumed that the system dynamics are described by a 
generic set of differential equations with unknown parameters. In this case the 
parameters in the dynamical model may be adjusted to ensure that the model 
reproduces available empirical data through a process of system identification. 
Alternatively the controller may be designed in such a way that it can cope 
with the uncertainty in the parameter values. 

4. Once the dynamical model of the plant is available, one develops (usually by 
analytical means) a control law that will modify the plant’s behavior to be in 
line with the original requirements. 

5. Finally the controllers are implemented and tested to ensure that the overall 
requirements are met. 

Much of control theory has been concerned with step 4, i.e. the synthesis of control 
laws. Ideally synthesis methods should: 

• be generic (i.e. can be applied in a straightforward fashion to most systems), 

• provide the designer with a convenient way to adjust the system’s response, and 

• guarantee stable system behavior. 

Many design methods have been developed and these mostly meet the objectives 
above, depending on the control objectives and also the class of system that is to be 
controlled. We can make the following broad (but not exhaustive) classification: 

A. Class of system: 

— Linear time invariant systems: i.e. systems with dynamics of the form: 

x = Ax + Bu (1.1) 
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where x is the state vector and u is the set of control inputs. 

- Linear time varying systems: i.e. systems with dynamics of the form: 

x = A(t)x. + B(t)u (1.2) 

Here the matrices A(t), B(t) are known functions of time. 

- Nonlinear systems: these are typically systems of the form: 

x = f(x,u), or (1.3) 

x = f(x) + G(x)u (1.4) 

Within the class of nonlinear systems there exists a special class of systems, 
i.e. globally feedback linearizable systems. These systems are relatively easy 
to control — see Section 1.1.1 for more details. 

- Nonlinear time-varying systems: i.e. systems with dynamics of the form: 

x = f(x,u,f) (1.5) 


B. Control objective: 

- Stabilization: Here the only objective is to find a control law that ensures 
that the system is (globally) stable. 

— Optimization: Here the objective is to find some control input that mini- 
mizes a cost function which typically has the form: 

J = x T {t f )Hx{t f )+ / \x T Qx + u T Ru)dT (1.6) 

Jto 

(It is implicitly assumed that part of the control objective is to ensure that 
the system remains stable) 

- Adaptive Control: In this case it is typically assumed that the system 
dynamics model contains some parameters, say p, which are unknown to 
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the designer, for example: 


x = f(x,p,u) (1.7) 

The controller then adjusts itself to accommodate the unknown parameters 
p and still maintain overall stability. 

— Robust control: Any of the above cases becomes a robust control problem 
if the objective is to maintain stability and performance objectives, de- 
spite the fact that assumptions regarding the nominal plant model may be 
incorrect. 

- Stochastic control: If the system dynamics in any of the cases above are 
affected by a random process we have a stochastic control problem. An 
example would be the following linear time varying system that is driven 
by a random process w: 


x = A(t)x + B(t)u + Tw (1.8) 

Table 1.1 gives a (very) brief summary of the availability of synthesis methods for 
different control problems. We note that the control of general nonlinear systems that 
are not feedback linearizable remains a very difficult problem. 

Our investigations will mostly be concerned with developing a synthesis method for 
general nonlinear systems which are not feedback linearizable. Since the global feed- 
back linearization theory also illuminates the general nature of the problem we want 
to address, we will start by examining it. Following that, we will briefly discuss op- 
timal control as an approach to synthesize control laws for nonlinear systems and 
then finally in Section 1.3 provide an overview of the methods to be developed in this 
thesis. 
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Linear 

Time 

Invariant 

Linear 

Time 

Varying 

Non-linear 

(feedback 

linearizable) 

Nonlinear 

(general) 

Stabilization 

Well developed 
[26] 

Via opt. 
control [32] 

Via LTI 
theory 

Unresolved 
(possibly via 
opt control) 

Optimization 

Well developed 
[32] 


Via LTI 
theory 

Hard to 
apply. [15] 

Adaptive Control 

Well developed 
[19] 

Partial results 
[43] 

Ongoing 
[57], [58] 


Robust Control 

Well developed 
and ongoing 
[60], [16] 

Via opt. 
control 

Ongoing 

[55] 

Unresolved 
(possibly via 
opt. control) 
[2] 

Stochastic- 

Control 

Well developed 
[32] 

Well developed 
[32] 

Not fully 
developed 

Hard to 
apply. [61] 


Table 1.1: Availability of Synthesis Methods 


1.1 Feedback Linearization of Nonlinear Systems 


Global feedback linearization and normal forms for nonlinear systems are concepts 
that have recently been developed and provide much insight into the fundamental 
issues regarding the control of nonlinear systems. These concepts have their origins 
in work done by Krener [31] and Brockett [7] among others. Further developments 
in the theory were provided by Hirschorn [22], Jacubczyk and Respondek [25] as well 
as Hunt, Su and Meyer [23]. In discussing these ideas we will focus on single-input 
single-output systems, although the concepts can be generalized for multivariable 
systems (see e.g. Slotine and Li [59]). 


1.1.1 Global Feedback Linearization 


We first consider systems which are so-called globally feedback linearizable. Such sys- 
tems are relatively easy to control since, in the words of Brockett, they are T . . linear 
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systems in disguise . . . More precisely the system: 


x = f(x) + g(x)tt (1.9) 

is considered to be globally feedback linearizable if a combination of coordinate change: 

z = <^(x) (1-10) 

and feedback: 


u = q(x) + /3(x)^ 

applied to (1.9) results in the linear time-invariant system: 

z = Az + bu 

where: 


0 1 0 ... 0 0 


0 

0 0 1 ... 0 0 


0 

0 0 

b = 


0 1 



0 0 0 ... 0 0 


1 


for all x. Technically it is required that: 


(1.11) 


( 1 . 12 ) 


(1.13) 


• the transformation: 

z = v?(x) (1.14) 

be a diffeomorphism, i.e. <^(.) must be smooth, and its inverse v 3_1 (-) must 
exist, and also be smooth. 

• the functions f(x) and g(x) be smooth. 

For a system of the type given in equation (1.9), precise conditions exist for us to 
determine whether a system is globally feedback linearizable or not (see e.g. Slotine 
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and Li [59], or Isidori [24]). In order to state these conditions compactly we use the 
following standard notation and definitions (see e.g. Slotine and Li [59]). Let h{x) be 
a smooth scalar valued function, and f(x) and g(x) be smooth vector valued functions 
of the n-dimensional vector x. 

• The gradient of h w.r.t x-coordinates is: 


Vh = 


dh 

dxi 

dh 

9x2 


dh 

dx n 


• The Jacobian of f w.r.t to x-coordinates is: 


Vf = 


(V/i) r 

(V/ 2 ) r 


(V /.) 1 


• The Lie-derivative of h w.r.t. f is given by: 

L f h = (Vh) T f 

and higher order Lie-derivatives are defined as: 

L° t h = h 

Lfh = LfL'f~ l h 

• The Lie-bracket of f and g is given by: 

adfg = [f,g] 

= Vgf-Vfg 


(1.15) 


(1.16) 


(1.17) 


(1.18) 

(1.19) 
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( 1 . 20 ) 

( 1 . 21 ) 



and higher order Lie-brackets are found from: 

ad?g = g (1-22) 

adjg = f,adf _1 g (1.23) 

• A linearly independent set of vector fields {fi,f 2 , • • • ,f m } is said to be involutive 
if and only if there exist scalar functions Q !J fc(x) such that: 

m 

= 52 a ufc(x)ffc(x) V*,j (1.24) 

k—l 

With the definitions above, the necessary and sufficient conditions for feedback lin- 
earizability of (1.9) can be compactly stated. We have the following powerful theorem 
(see e.g. Slotine and Li [59] for details): 

Theorem 1.1.1 The nonlinear system (1.9) is globally feedback linearizable if and 
only if the following conditions hold for all x; 

• The vector fields 

{g,ad f g adT'g} (1.25) 

are linearly independent. 

• The vector fields 

{g,ad f g,...,adr 2 g} (1-26) 

are involutive. 

In proving the theorem above one also obtains a method to construct the coordi- 
nate transform and feedback law required to put (1.9) into the linear form of equa- 
tions (1.12) and (1.13). The construction, which we describe next, requires one to 
solve a potentially difficult set of partial differential equations. 
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Assuming that conditions (1.25) and (1.26) are satisfied, the process to feedback 
linearize the system will be: 


1. Solve the following set of p.d.e’s for .pi : 


V^fadfg = 0 z = 0, .... n — 2 

V^[ad^ _1 g ± 0 


2. Apply the coordinate transform: 


z = ¥>(x) = 




L 


n— 1 

f 


r 5 ! 


which will result in the following dynamics: 


2l 


Z-2 

k 

— 

^3 

. * n . 


q(z) + /3(z)u 


where: 


<*Mx)) = Lftpi 
/3Mx)) = ig^rVi 


(1.27) 

(1.28) 


(1.29) 


(1.30) 


(1.31) 

(1.32) 


We shall refer to (1.30) as the cascade canonical form 1 because the system is 
now made up of a string of integrators with the highest derivative driven by a 
nonlinear function of the state as well as the control input. 

'The author is not aware of a standard term used for this canonical form for nonlinear systems 
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3. Apply the feedback: 


“ = ( -" Wx » + ^ 


(1.33) 


Once we have applied the feedback and transformations above we can apply the large 
body of control theory for linear systems. However the process outlined above is 
easier stated than applied. Difficulties arise on the following counts: 


1. Establishing that a system is globally feedback linearizable (so that we know 
that the p.d.e’s (1.28) have a solution) can require significant effort. For in- 
stance, if we cannot analytically check the conditions of Theorem 1.1.1, we 
have to test these conditions numerically at “every” point in the state space. 

2. We have to solve the partial differential equations (1.28) which may be a far 
from trivial problem to do either analytically or numerically. 

3. Once we have applied “inner loop” control (1.33), we expect to apply further 
state feedback control to get the L.T.I. system to behave appropriately, e.g.: 


= k r z 


^1 ? ^2 1 , kn 


L f ipi 

i 


(1.34) 


(1.35) 


For both the case where we can find z = <^(x) analytically as well as the case 
where we find ifi numerically, the feedback law might be quite sensitive to 
errors in x. Since the “outer loop” feedback. v> = k T z, involves derivatives of 
yd we would expect any problems to be exacerbated if we have to find ^ by 
numerically solving the p.d.e.’s (1.28). 
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1.1.2 Partial Feedback Linearization and Normal Forms 

In the previous section we outlined the situation where a system can be globally feed- 
back linearized. As we might expect, there exists a large class of systems which cannot 
be feedback linearized, and thus remain difficult to control. To better understand the 
issues involved the concept of partial feedback linearization will be examined next. 
We shall see that if f(x) and g(x) are smooth enough, systems of the form (1.9) can 
be (locally) partially feedback linearized and put in a normal form. 

Consider controllable systems of the form: 


x = f(x) + g(x)u (1.36) 

y = h(x) (1.37) 


where y = h(x) is an arbitrary smooth function of the system state. Usually y will be 
some output variable that we want to control. In its simplest form partial feedback 
linearization is achieved as follows: We successively differentiate y until the control 
input appears in the expression for say the rth derivative of y, viz: 


dr 


d /d r l y 
dx i d r -1 


(f(x) + g(x)u) 


= Lfh + L % L\~ x h 
= ct(x) + l3(x)u 


The state equations for the system can then be written as: 


y i 


y 2 

y 2 


y 3 

Vr 


a(x) + /3(x)u 



1 

^■T 

U/ 

1 


(1.38) 

(1.39) 

(1.40) 


(1.41) 
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If r = n, where n = is the dimension of the state vector, and we use the control: 

« = -q7~\ (~ Q (x) + v) (1.42) 

P( x ) 

we have globally feedback linearized the system. The process for global feedback 
linearization can thus be thought of as a method to appropriately construct y = h(x), 
so that r = n, if it is at all possible. Furthermore, the conditions of Theorem 1.1.1 
can then be thought of as tests to determine whether a suitable h(x) can be found. 

In equation (1.41) we see that ^(x,u) depends on both x and u. It can be shown 
that using an appropriate transformation (which is also constructed by solving partial 
differential equations) the system can be put in the following normal form: 

2/2 
2/3 

<*(y,v) + @(y,v)u 

w(y,i?) 

This normal form enables us to better understand what we can achieve in controlling 
the system. For instance, if /?(y,tj) 0 Vx, we can fully determine the input-output 

relation between u and y . However, since y acts as an input to the 77 states, we 
see that y following a certain trajectory forces a specific behavior for the internal 
dynamics: 

^ = w(y,rj) (1.44) 

Clearly it is possible that the states r /, which are not observable via h(x ), can become 
unstable, depending on the y trajectory. 

An important special case occurs when y = 0 . The dynamics of (1.43) subject to the 
constraint y = 0 are called zero-dynamics [59] of the system. In terms of the normal 


( 1 . 43 ) 
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form we get: 


y 

v 


0 

w (0, 17) 


(1.45) 


If the zero-dynamics of the system is asymptotically stable the nonlinear system is said 
to be asymptotically minimum phase. It can be shown that this definition captures 
our standard understanding of minimum-phase systems for the case where w'e are 
dealing with linear systems (see Slotine and Li [59]). 


At this point it becomes clear what the challenge is regarding the control of nonlinear 
systems. Typically the problem will be to obtain certain desirable behavior from our 
output variable y while still maintaining stability of the internal dynamics. 


1.2 Optimal Control 

Another approach to synthesizing controllers for nonlinear systems is to find the 
control input through solving an optimal control problem. In principle this approach 
meets all our objectives for an ideal synthesis method. For example if our objective 
is to stabilize a nonlinear system of the form: 

x = f (x, u) with equilibrium point: ( 1 .46) 

0 = f(0,0) ( 1 . 47 ) 

We then pose an optimization problem w r hich has the objective to minimize the fol- 
lowing cost function: 

roc 

J= {x T Qx + u T Ru)dr (1-48) 

j to 

where we assume that both Q and R are positive definite matrices. If the system (1.46) 
is controllable [63] and we apply the control input that minimizes (1.48) the system 
will be stable as the following standard argument shows: 

Since (1.46) is controllable there exists some control input that will drive the system 
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to the origin from any initial condition. If we apply this control input to the system 
we will incur a certain finite cost, say J c (x o). Furthermore if we apply the optimal 
control input to the system there will be a certain associated cost, say J(x 0 ) which 
has to be less than J c and thus finite. Now note that since Q is positive definite the 
cost incurred for an unstable response would be infinite. Hence the optimal input 
cannot cause unstable responses. 

A further attractive aspect of the optimal control approach is that it provides us with 
a method to conveniently adjust the system responses by adjusting the cost function. 

The only undesirable feature of the optimal control approach is that it is compu- 
tationally very demanding. If the goal is to obtain the optimal control as a state 
feedback law we have to solve a nonlinear partial differential equation, the Hamilton- 
Jacobi-Bellman equation , for the whole region of the state space where the system is 
to operate (see Section 3.2.1 for further discussion). Alternatively we have to find the 
optimal control input by numerically solving a two point boundary value problem [15]. 


1.3 Overview 

As indicated earlier we will be mostly concerned with developing a controller synthesis 
method for nonlinear systems of the form: 

x = f(x) + G(x)u (1-49) 

We will exploit the fact (see Chapter 2) that provided f(x) meets certain smoothness 
requirements, equation (1.49) can be written in the following quasilinear form: 

x = A(x)x + B(x)u (1.50) 

In Chapter 3 we will investigate a design methodology where the control input is com- 
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puted by continuously solving a matrix Ricatti equation based on the instantaneous 
properties of A(x), B(x). We shall refer to this as a control law with zero- look-ahead. 
Related work appears in the thesis by Shamma [53] who investigates gain scheduling 
in the context of the LQG\LTR methodology [60], and takes a different perspective 
with regard to the stability issue. 

In Chapter 4 we develop control laws that take into account the “future behavior" 
of ,4(x),£?(x). We refer to these as control laws with look-ahead. The methods we 
investigate will all have underlying optimal control problems that have to be solved. 
By exploiting the quasilinear form we will obtain approximate solutions to these 
optimal control problems. 

Finally in Chapter 5 we will apply the controller synthesis methods we developed to 
control a complex nonlinear system, i.e. a flexible robotic manipulator similar to the 
space shuttle Remote Manipulator System [20]. 


25 



Chapter 2 


Quasilinear Description 


In this chapter we will examine some of the basic issues in dealing with nonlinear 
dynamical systems of the form: 


x = y4(x)x + B(x)u 

We shall refer to such systems as quasilinear systems. 


( 2 . 1 ) 


2.1 Generality of the Quasilinear Form 

The following lemma from Vidyasagar [63] shows that under mild differentiability 
conditions, systems of the form: 


x = f (x) + G(x)u (2.2) 

can always be transformed to the quasilinear form of equation (2.1). The lemma is 
more general than we will need since it caters for the case where f(0) ^ 0, whereas 
we will generally assume f(0) = 0. 
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Lemma 2.1.1 Suppose f : R n — > R n is continuously differentiable. Then there 
exists a continuous function A : R n — » R nxn such that: 

f (x) = f(0) + .4(x)x. Vx € R n (2.3) 


Proof: 

Fix x and consider f(Ax) as a function of the scalar parameter A. Then 


f(x) = m+l ^f(Ax)dA 

= f(0) + [jfv x f(Ax)dA 


so that 


A(x)= f l V x f(Ax) dX 
Jo 


(2.4) 

(2.5) 


(2.6) 


Q.E.D. 


The above formula gives us a specific quasilinear form for any C 1 function f(x). 
Equation (2.6) is useful in cases where the quasilinear form may not be immediately 
apparent. For example let: 


f(x) 


sin(xi + x 2 ) 
sin(x 2 ) 


2.7 


Equation (2.6) gives: 


*(X) 


sin(ri +J2 ) sin(ri +-T2) 
(^1+^2) (xi+x 2 ) 

0 sin( -r 2 ) 

x 2 


( 2 . 8 ) 


Note that, in general, the quasilinear form for a given function f(x) is not unique. 
In fact, we see that we can add any vector that is orthogonal to x, to any row of 
,4(x), and still get the same value for A(x)x. For the example above we see that an 
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alternative quasilinear form for f(x) is: 

£1 

x 2 _ 

This non-uniqueness can be exploited as we shall see in Section 2.3.2. 


f(x) = 


sm(x,+x 2 ) _ k sin (*,+*,) 

Xi+X 2 Xi+X 2 1 


— k 2 X2 


sin(x2) 

^2 


+ ^2X1 


(2.9) 


2.2 Quasilinear Form for Mechanical Systems 


In this section we will show that the equations of motion of a class of mechanical 
systems which operate in the absence of gravity can be put into quasilinear form 
without solving the potentially difficult integrals of equation (2.6). (See also Ap- 
pendix A which describes a set of computer routines that utilize this approach to 
automatically derive quasilinear equations of motion for flexible manipulators.) 


We will use Lagrange’s equations. They are: 


d.dL dL 

7t { ai ] ~ % = Q " 


where: 

n = the number of degrees of freedom in the system 
qi = ith generalized coordinate 
Q, = ith generalized force 
L = T — U 
T = kinetic energy 
U = potential energy 


( 2 . 10 ) 


If we can express the kinetic and potential energy of the system in the form: 


T = iq r W(q)q (2.11) 

U = t q r A'q (2.12) 
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where: 


q r = [?i, ? 2 ; • • • , q n \ is a vector of generalized coordinates 
H( q) = a generalized inertia matrix 
K = a generalized stiffness matrix 
then Lagrange's equations can be written in vector format as: 


H q + H q - V q T + Aq = Q 


(2.13) 


where: Q = \Q\, Q 2 , • • • . Qn] T is a vector of generalized forces, and: 


r\ 

v q r = —t 

dq 

'dT_ dT_ dT_ 
~ dq\ dq 2 ' dq n 


Now: 


v q r = 


q T l^(//( q )q) 


d 


iT 


so that the equations of motion can be expressed as: 


H q + 




lT\ 


H{ q)q 


q + A"q = Q 


(2.14) 


(2.15) 


(2.16) 


(2.17) 


or 

#(q)q + C(q,q)q + Aq = Q (2.18) 

To get a state space description we furthermore assume that the system is controlled 
via the generalized forces and that they can be expressed as: 


Q = M(q,q)u (2.19) 

where the components of u are the control inputs to the system. We then get the 


29 



desired state space form: 


q 

q 


o i 

-H(q)-'K -//(q)- 1 C(q,q) 


q o 

+ U (2.20) 

qj [//(q)' 1 M(q,q)_ 


Note that the key to getting to the quasilinear form is to express the kinetic and 
potential energy as in equations (2.11) and (2.12). In the following example we shall 
use a double pendulum system to show how we can find the desired form for the 
kinetic and potential energy. We shall see that it is generally relatively easy to find 
the desired form for the kinetic energy. On the other hand, the expression for potential 
energy will only be in the desired form if we are dealing with the equivalent of “linear 
springs’ 1 . This will typically be the case for flexible mechanical space based systems. 
If we are dealing with nonlinear “springs”, e.g. gravity, we generally will have to use 
equation (2.6). 


Example 2*2*1 

Figure (2-1) shows a double pendulum with rigid massless links of lengths l\ and / 2 
respectively. The links also have tip masses mi and m 2 and rotational springs K \ 
and K 2 attached. For this system the generalized coordinates are 9 X and 0 2 , viz: 



( 2 . 21 ) 


We first find an expression for the kinetic energy. Let Ri and R 2 be position vectors 
to each of the masses. Then kinetic energy of the system is given by: 


mi /dRi dRi\ m 2 /dR 2 dR 2 \ 

2 \ dt d£ ) + 2 \ dt d t ) 


with: 


Ri = 


/ 1 cos(^i) 
li sin($i) 


, r 2 = 


l\ COs($i ) + / 2 COs({?i + 62) 
/ 2 sin(#i ) -j- / 2 sin(^i -(- @ 2 ) 


For systems not made up of point masses we can simply express the kinetic energy 
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Figure 2-1: Double Pendulum 


as an integral over infinitesimal mass elements viz: 


T = 


r dK 

Jv oi d£ d t 


dR dR 


dm 


For the double pendulum system we get: 


dR, 



Ox 

d t 

9 D 9 p 

d&i n '» d0 2 rLj 


0 2 


= J,(q)q 


where: 


^i(q) 


—l\ sin(0i) 0 

l\ 008 ( 0 !) 0 


(2.24) 


(2.25) 


(2.26) 


(2.27) 
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I" — /i sin(0i) - /'2 sin(^i + 0 2 ) ~h sin(£h + 0 2 ) 


M q) = = (2-28) 

lx cos(^i) + / 2 cos($i + 0 2 ) h cos(0i + 0 2 ) 

Now we can express the kinetic energy in the desired form: 

T = ^q+jq^q (2.29) 

= iq T [rmJfjj + m 2 jJ J 2 ] q (2.30) 

= ^q T //(q)q (2-31) 


Next we find an expression for the potential energy. In the absence of gravity the 
potential energy of the system is given by: 


= \(K x el + K 2 6 2 ^ 

(2.32) 

1 T [ K x 0 

(2.33) 

= q 

2 0 k 2 


1 T 

m jq T A-q 

(2.34) 


which is of the desired form viz. equation (2.12). We achieved this form because 
we had linear springs. If we include the effects of gravity we will not be able to 
conveniently express the potential energy in the form of equation (2.12) - gravity 
does not act like a linear spring. Nevertheless we can still obtain a quasilinear form 
for our dynamics when gravity is present by e.g. using (2.6) as we shall now show. 

The potential energy due to gravity is given by: 

U gravity = ~m x gli cos(6 l 1 ) - m 2 g (/1 cos(^) + l 2 cos(0 x +0 2 )) (2.35) 

This gives rise to the following additional term in equation (2.13): 

m\gl\ sin(0 t ) — m 2 g {—l\ sin($i) — l 2 sin(0 x + 0 2 )) 
v q V gravity ~ (2.36) 

m 2 gl 2 sin(0! + 0 2 ) 
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In this case we can express the gravity term in quasilinear form as follows: 


y r 

v q ^'gravity 


(mi/i+m 2 /i)ff sin(fli ) , m 2 g/ 2 sin(fli ±$ 2 ) 

0i (6>i+0 2 ) 


m 2 g/2 sin(^i +fl 2 ) 

(0i+0 2 ) 


G 


3 r av 


q)q 


m 2 gl 2 sin(gi+g 2 ) 
( Q \ +#2 ) 

m 2 gl2 s»n(#i +9 2 ) 

(#i +# 2 ) 



’ 6i 


1 

CN 

<3S 

i 


(2.38) 


We then get a result similar to equation (2.17): 


H q + 



_d_ 


1 T\ 


m q)q 


q+ [A + G 5rav (q)] q — 0 


(2.39) 


The rest of the process follows as outlined above. The resulting quasilinear equations 
of motion are given in Appendix B. 


Remarks: 


• The quasilinear form of equation (2.20) is already in a partially feedback lin- 
earized form, although not the normal form of (1.43). However, we are using 
"natural” coordinates, i.e. the variables of the equations of motion have physi- 
cal significance, whereas the variables rj in the normal form of (1.43) usually do 
not have physical significance. 

• If the number of actuators is equal to the number of degrees of freedom of the 
system, the quasilinear equations will be in the desirable cascade canonical form 
of equation (1.30). This will, for instance, be the case if we are dealing with 
robots with rigid links with actuators at each joint. 


2.3 Stability Analysis Using the Properties of A(x) 

In this section we will examine some stability tests for the ordinary differential equa- 
tion (o.d.e.): 

x = A(x)x (2.40) 
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that are based only on the properties of the matrix A(x). Because the o.d.e. (2.40) 
represents a very large class of nonlinear systems, viz. any o.d.e. x = f(x) with 
f(.) £ C 1 , we do not expect to find stability tests that are easy to apply and at the 
same time provide both necessary and sufficient conditions for stability. We note 
that the tests we discuss will also be useful for analyzing the stability of linear time 
varying systems where: 

x = A(t)x (2.41) 

2.3.1 Stability Analysis Based on Eigenstructure of A(x) 

One is tempted to predict the stability of (2.40) based on the eigenvalues of A(x). We 
will show that the eigenvalues do not provide us with enough information to deduce 
stability and that we have to take into account the rate of change of the eigenstructure 
of A(x). The results we present are similar in spirit to other results regarding the 
stability of slowly varying systems (see e.g. [43], [63]). 

We have 

x = A(x)x (2.42) 

Assume that all the eigenvalues of A(x) are distinct, real and strictly negative. Then: 

A(x) = r(x)A(x)7’~ 1 (x) (2.43) 


where 


T(x) = a matrix containing the eigenvectors of A(x) 

A(x) = diag(Aj(x), A 2 (x), . . . , A n (x)) is a diagonal matrix containing the 
eigenvalues of A(x) with A;(x) < — t < 0 Vx, i = 1 . . . n. 
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We can then make the coordinate change: 


z = T 1 (x)x (2.44) 

which results in the following dynamics for the transformed system: 

z = T _1 (x)x+ ^ (r _1 (x))x (2.45) 

= T~ 1 A(x) (rr _1 ) x — T~ 1 TT~ 1 x (2.46) 

= A(z)z - (T _1 (z)jT(z)) z (2.47) 


where for notational convenience we have not shown the state dependence of T in 
equation (2.46). To determine the stability of the system we use the following simple 
observation: 

Fact 2.3.1 The system of o.d.e.’s: 


z = A(z)z + M(z)z 

with: 

Ai(z) 0 ... 0 

0 A 2 (z) ... 0 

A(z) = 

0 0 ... A n (z) 

will be asymptotically stable provided that: 

z r A(z)z + z r M(z)z <0 Vz ^ 0 


Proof: 

Use the candidate Lyapunov function: 


U(z) 



Z 


(2.48) 


(2.49) 


(2.50) 


(2.51) 
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which is positive definite and radially unbounded. Then: 


V = z T z (2.52) 

= z t A(z)z + z t M(z)z (2.53) 

and global asymptotic stability follows since V < 0 Vz / 0. 

Q.E.D. 


We can apply this fact to our case by noting that: 

z t A(z)z < -ez T z (2.54) 

and 

z t T-'(z)T{z)z < ||z||, ||r- , (*)f(*)*|| 1 
< ||r“‘(z)T(z)|| 2 z T z 

so that if; 

||r' 1 (z)T(z)|| 2 < u < e (2.57) 

we have: 

z r A(z)z + z t T~ 1 Tz < (-e + v) z T z < 0 (2.58) 

That is if the rate of change of the eigenstructure is small enough, viz. T « 0, we 
can infer stability pointwise from the eigenvalues of A(x). For example, consider the 
following seemingly complex second order nonlinear system: 

x' x = —(15x2 + 1260 xjX 2 + 1890xjX2 + 1890xjX2 

+ 1260 XjX 2 + 540 xjX 2 + 135 xi x® 

+ 135x®x 2 + 540 x\x\ + 51 x 2 + 54 Xi + 15 x^ 

—6 x\ — 15 xx x 2 — 9 x 2 — 27 XjX 2 — 45 x^x 2 


(2.55) 

(2.56) 
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-33 x x x 8 + 461 x\ x 2 + 550 xi x\ + 432 x*x 2 
+ 1026 xiX 2 + I2O6X1X2 + 702 x! X2 - 6x1 
—9 x 2 + 130 Xj + 219 Xj + 72 Xj + 162 x\ + 15 x\ 

+ 15x2 + 105 xjX 2 + 315xiX2 + 525xjX 2 + 525 XjX* + 
315 xjX 2 + 105 xj x 8 )/ ( 1 + xj + 2xi x 2 + x\) 
x 2 = ( 10 x^ + 840 x?x 8 + 1260 x\x\ + 1260 x\x\ 

+840 XjX® + 360 xjX 2 + 90 x x x 8 + 90xjX 2 
+360x(x 2 + 31 x 2 + 34x! + 10xj — 4x 2 — 10xi x 2 
— 6 x 2 — 18xiX 2 — 30xjX 2 — 22 x t x\ + 305x 2 x 2 
+364 Xi x^ + 288 x|x 2 + 684 XjX^ + 804 x 2 x^ 

+468 xi X 2 — 4 Xj — 6 x^ + 86 Xj + 145 x 3 2 + 48 Xj 
+ 108 x* + 10 x] + 10 x' + 70 x®x 2 + 210 x\x\ + 350 x?x 8 
+350x 3 x 2 + 210xjX 2 + 70x x x 2 )/(l + x 2 + 2x! x 2 + x'l) 

These equations can be written in quasilinear form as: 

x = A(x)x 


where 


A(l,l) = -{54 + 58 (x! + x 2 ) 2 + 18 (-2xj -3 x 2 ) 2 
+ 18 (—2 xi — 3 x 2 ) 2 (xi + x 2 ) 2 
6 xi — 9 x 2 + 3 ( — 2 xi — 3 x 2 ) (xi + x 2 ) 2 
+ 15 (xi + x 2 ) 6 + 15 (xi + x 2 ) 8 j / {l + (xi + x 2 ) 2 j 
21(1,2) = -{51 +57 (xj + x 2 ) 2 + 18 (-2xj -3 x 2 ) 2 
+ 18 (— 2 Xi — 3 x 2 ) (xi + x 2 ) 2 
6 x 1 9 x 2 + 3 ( — 2 xi — 3 x 2 ) (xi + x 2 ) 


(2.59) 


(2.60) 


(2.61) 


(2.62) 
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+ 15 (#i + x 2 ) 6 + 15 (xi + ^2) 8 | / |l + (*Ti + x 2 ) 2 | (2.63) 

A(2A) = {34 + 38 (xi + x 2 ) 2 + 12 (— 2xi — 3x 2 ) 2 

+ 12 ( — 2 x 1 — 3 x 2 ) (xi + x 2 ) “ 4 x 1 — 6 x 2 
+2 (— 2 Xi — 3 x 2 ) (xi + x 2 ) 2 + 10 (xi + X 2) 6 
+ 10 (xi + X 2 ) 8 } / jl + (xi + x 2 ) 2 } (2.64) 

A( 2 1 2) = {31 + 37 (xi + X 2 ) + 12 (— -2 x\ — 3 X 2) 2 
+ 12 (—2xi — 3 x 2 ) 2 (xj + x 2 ) 2 — 4xi 
-6x 2 + 2 (-2xi - 3 x 2 ) (xi + x 2 ) 2 + 10 (xi + x 2 ) 6 
+ 10 ( x 1 + x 2 ) { / { 1 + ( x 1 + x 2 ) { (2.65) 


For A(x) above, an eigenvector matrix is: 


-1 -3/2 
1 1 


( 2 . 66 ) 


Since eigenvectors are only uniquely determined up to a scale factor we use the fol- 
lowing (eigenvector) transformation matrix: 



where we assume that a and j3 are constants. For these eigenvectors the eigenvalues 
of A(x) are given by: 


Ax(x) 

A 2 (x) 


3 + x\ + 2 £2 + C9 fi8'l 

1 + xj + 2 a?! x 2 4- x\ 

(—5 x 2 — 30 x\ x\ — 75 x\x\ — 100 — 54 x\ — 75 x\x\ 

—30 x\x 2 + 3 x 2 — 72 xx x 2 — 24 x\ — 20 + 2 x\ — 5 + (2.69) 
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Then upon applying the state transform: 


z = r~ ! x 


(2.70) 


and noting that T = 0 we get: 

z = A(z)z (2.71) 

with: 


A(z) = 


-(12 + ^f)/(4 + /^f) 
0 


0 


-20 + a z i—6 q 2 ^j 




We see that Aj(z) < 0, A 2 (z) < 0 and also T = 0 and hence conclude that the system 
is stable. 


2.3.2 Stability Analysis via Matrix Measure Theorems 

Another approach to inferring stability of (2.40) from the properties of .4(x) is to use 
the matrix measure or logarithmic derivative of A(x) which is given by: 


Definition 2.3.1 Let A be an n x n matrix and ||.|| t - be an induced norm on R nXn . 
The corresponding matrix measure p{A) of the matrix A is given by: 


MA) 


lim 

o+ 


||7 + eA.||t - 1 

t 


(2.73) 


The matrix measure was originally used to determine stability and error bounds for 
numerical integration methods [13] and later used by Coppel [12] to obtain solution 
bounds for linear time varying o.d.e.’s. Desoer and Haneda [14] also used matrix 
measures to bound solutions of o.d.e.’s, bound accumulated truncation errors for 
numerical integration schemes and determine convergence regions for the Newton- 
Raphson solutions of nonlinear equations. Further related work appears in [51]. [52] 
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and [44]. 

The matrix measure is defined so as to easily bound the behavior of ||x(t)|| for o.d.e.'s 
of the form: 

x = A(x,t)x (2.74) 

The following theorem shows how this is done. 

Theorem 2.3.1 [63], [14] Consider the o.d.e.: 

x = A(x,t)x, x € R n (2.75) 

Let ||.|| be a norm on R n , and ||.||, and p(.) be the corresponding induced norm and 
matrix measure on R nxn respectively. Then: 

IjIMOII < M4(x,())||x(()|| (2.76) 

Suppose furthermore that there exist continuous functions a(t),j3(t) such that: 

p(A(x,t)) < a(t), /3(t) < fx( — A(x,t)) V£>(LVx£i? n (2.77) 

Then: 

M<o)l|exp (f -P{r)d r) < ||x(<)|| < ||x(f 0 )||exp ^ a(r)dr) (2.78) 

Proof: 

We first show (2.76). We have from (2.75): 


x(t + St) = x(t) + A(x,t)x(t)6t + o(6t) (2.79) 

= [I + 6tA(x, t)]x(t) + o(St) (2.80) 

=H|x(t + £f)|| ^ || [I + St A(x,t)]x(t)|| + ||o(<5t)|| (2.81) 

< ||/ + «M(x,f)yx(«)|| + ||o(ft)ll ( 2 . 82 ) 
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|x(< + «f)ll-||x(OII < {||/ + 5 M(x,0||,--1}||x(0IH-|K^)II (2.83) 

| -' in , 


j«»> * ,& r‘T 


Equation (2.78) then follows by multiplying both sides of the inequality: 


d+ 

dt 


|x( # ) || < lim 


\\I + ^Afx 


st~ o+ ( 

< a(0l|x(0ll 




(Oil 


(2.85) 

( 2 . 86 ) 


with the integrating factor: 


exp 


J to o(r)dr^ 


( 2.81 


The matrix measures associated with the 1,2 and oc norms are given by: 


j Q.E.D . 


||x||oo = max(xi) 

t 

l|x||i = W 

i 

Mb =(£>,■ 


2=1 

" A ,/! 


<1=1 


// 00 (A) = max Ui, + V] |a,j| 

\ 

Vi(A) = max o^A^KI 

3 \ 

// 2 (A) = \ max (A t + a) /2 


( 2 . 88 ) 

(2.89) 

(2.90) 


The usefulness of matrix measures becomes apparent by examining equations (2.77) 
and (2.78) and noting that //(A) can have negative values and thus be used e.g. to 
show that ||x|| — ► 0. The following example shows how this can be done: 


Example 2.3.1 Consider the following set of o.d.e.’s: 


i\ = — (3 + X!X 2 ) 2 Xi + cos(x x ) sin(x 2 )x 2 (2.91) 

x 2 = xi sin(xj) — x 2 Xj — 2 x 2 (2.92) 
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These can be put in quasilinear form as: 


i\ 


— (3 -f X 1 X 2) 2 cos(xi) sin(x- 2 ) 


Xi 

x 2 


sin(xj) —2 — X 2 


x 2 


For A(x) above we have: 

Pi(T(x)) < -1 

/ioo(A(x)) < -1 


(2.93) 


(2.94) 

(2.95) 


Note that the quasilinear form is non-unique (see Section 2.1), and that the realization 
above is much more useful than the alternative: 


x x 


— (3 + xix 2 ) 2 cos(xi) sin(x 2 ) 


X\ 

X‘2 


sin(xj) — X 1 X 2 —2 


*2 


(2.96) 


where now the —X\X 2 term in A( 2, 1) can become unbounded and thus make it im- 
possible to establish the row or column dominance of the diagonal terms as required 
by equations (2.88) and (2.89) 


Remarks: 


• From equation (2.90) we see that the system (2.75) will be stable if A(x) is 
symmetric and negative definite. 

• Equations (2.88) and (2.89) show that in order to have ^(A(x)) < 0 (which 
we need for stability, see equation (2.78)) the diagonal terms of .4(x) have to 
be negative and dominate the off-diagonal terms. Unfortunately, when dealing 
with physical systems we often naturally obtain dynamical equations that do 
not result in A(x) being diagonally dominant as we see in the following example 
and is also apparent in Section 2.2. 
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Figure 2-2: Nonlinear Spring Mass Damper 

Example 2.3.2 Consider the spring mass damper system shown in Figure (2-2). The 
spring and damper are both nonlinear with stiffness g(z) and damping characteristic 
f(x,x) respectively. The dynamics are governed by: 


x + f(x,x)x + g(x)x = 0 


(2.97) 


These equations are then conveniently put into the following quasilinear form: 


x 

x 


0 1 

-g(x) 


X 

X 


(2.98) 


Here we see that for the 1 and oc-norm cases we are stymied by the fact that ,4(1, 1] 
0. Using the ‘2-norm does not help either because we see that: 


det (A/ - (A t + 24)) = A 2 + 2/A - (g - l) 2 


(2.99) 


which has roots at: 


c-4 +i ^T 


at least one of which will be non-negative. Hence: 


( 2 . 100 ) 


/i(.4(x)) > 0, Vx 


(2.101; 
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The challenge seems to be to find a transformation z = T(x)x such that in the 
transformed system the matrix A(x) will have negative entries on the diagonal that 
dominate the off-diagonal terms. If A(x) is in phase variable form and its eigenvalues 
change relatively slowly we can use a Vandermonde matrix to effect the transform. 
For instance, given: 


A(x) 


0 1 0 ... 0 

0 0 1 ... 0 

0 0 0 0 
0 0 0 ... 1 

a 0 (x) ai(x) a 2 (x) ... a„_ x (x) 


( 2 . 102 ) 


with eigenvalues A x (x), A 2 (x), . . . A n (x). Then using the coordinate transform: 




z = 

T~\x)x 


(2.103) 

where: 

1 

1 

1 

1 



Aj(x) 

A 2 (x) 

A 3 (x) 

A„(x) 


r(x) = 

Mx ) 2 

A 2 (x) 2 

A 3 (x) 2 . 

... A„(x) 2 

(2.104) 


_ A^x )- 1 


A 3 (x) n_1 . 

3 

£ * 

1 

i 


we get: 








z = A(z)z - T~ x Tz 


(2.105) 

so that if the eigenvalues change relatively slowly we 

have T^O and the transformed 


system will be diagonally dominant. 
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2.4 Summary 


In this chapter we established the following: 

• A broad class of dynamical systems of practical importance can be represented 
in quasilinear form. 

• The dynamics of space-based flexible mechanical systems can easily be expressed 
in quasilinear form. 

We also showed some sufficient conditions for the stability of: 

x = A(x)x (2.106) 

based on the properties of A(x). In particular we noted that: 

• If A(x) had negative real eigenvalues bounded away from the imaginary axis 
and the rate of change of the eigenvectors was small enough the system would 
be stable. 

• If any matrix measure of A(x), say //(A(x)), is “essentially negative”, i.e. 

lim^ooexp (/ 0 * A(x(r)))dr^ = 0 the system will be stable. Such cases would 

for example occur when: 

— A(x) had negative diagonal entries which dominated the off-diagonal terms. 

— A(x) was symmetric and negative definite. 

We should emphasize that these tests provide only sufficient conditions for stability 
based on the properties of A(x) alone and that systems may well be stable even when 
these tests fail. 

In the next chapter we will show how we exploit the quasilinear form to develop con- 
trollers for some nonlinear systems, and also investigate further the stability question 
using Zubov’s method. 


45 



Chapter 3 


Quasilinear Controllers with Zero 
Look-ahead 

3.1 Introduction 

As indicated in Chapter 1, we would like to have a controller synthesis method that 

• is generic (i.e. can be applied in a straightforward fashion to most systems) 

• guarantees global stability 

• provides the designer with a convenient way to adjust the system’s response. 

We also noted that there are few methods that meet these requirements. An attrac- 
tive method was that of global feedback linearization but it could only be applied to 
a limited class of systems. Furthermore we saw that even when systems were glob- 
ally feedback linearizable, the process of actually obtaining the feedback laws could 
still be difficult. Another approach was to synthesize control laws by solving an opti- 
mization problem. We noted that in principle the optimization approach satisfied our 
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requirements for an ideal synthesis method but could be computationally demanding. 
In this chapter we will consider a method for control system design that is related to 
the optimal control approach but requires less computation. We will refer to it as the 
continuous Ricatti design method (see also [53] for related work). 

The method we introduce exploits the quasilinear form discussed in Chapter 2. It 
has been successfully applied to nonlinear systems in simulation studies and partially 
meets our requirements for a control synthesis method, i.e. it is generic in nature 
and provides the designer with design parameters to conveniently influence system 
responses. However only local stability can be guaranteed. The stability domain 
for the closed loop system will have to be determined “after the fact” (see also Sec- 
tion 3.4). 


3.2 The Continuous Ricatti Design Method 

We develop the continuous Ricatti design method using the Hamilton- Jacobi- Bellman 
(H-J-B) theory [28] for optimal control as a starting point. By utilizing the quasilin- 
ear form we will find a candidate solution to the H-J-B-equation and an associated 
feedback law. We will use this associated feedback law to actually control nonlinear 
systems despite the fact that it may not be optimal (non-optimal because it is found 
from a candidate solution to the H-J-B equation). 

3.2.1 The H-J-B-Equation and Quasilinear Dynamics 

H-J-B-theory deals with the following nonlinear optimal control problem: 

Find the control input u that minimizes the cost function: 

J = / / £(x(r),u(r),r)dr (3.1) 

J to 
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where: 


£(x(r),u(r),r) > 0 
and the system dynamics are given by: 

x = f(x,u,f) 
x(t 0 ) = x 0 


(3.2) 


(3.3) 


(3.4) 


To solve the optimization problem above H-J-B theory uses the imbedding principle 
and considers the problem above as a specific case of a larger set of problems. The 
larger set of problems is constructed by allowing the initial condition to take any 
admissible value (instead of only x 0 ) and the initial time to take any value between 
t 0 and tj. Thus the cost function for the larger set of problems becomes a function 
of the variable initial state, say x, and variable initial time, say t , viz: 

J(x,<) = j ' £(x(r),u(r),T)dr (3.5) 


Assuming that a solution to the optimization problem exists it can be shown (see 
e.g. [28]) that the optimal control for the original problem (3.1)-(3.4) is given by: 

u opt = argmin j(£(x,u,r)) + (f(x,u,f))j (3.6) 

where J(x,t) is found by solving the Hamilton- Jacobi- Bellman partial differential 
equation (H-J-B-equation): 


<9J(x, t) 
dt 


— min 

u 


(£(x, u, r )) + 


dj(x , t ) 

dx 


(f(x,u ,<)) 


1 


(3.7) 


subject to the boundary conditions: 


J(0,t) = 0 


(3.8) 
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J(x, t) > 0 Vx^O 


(3.9) 


It can furthermore be shown that these conditions provide both necessary and suffi- 
cient conditions for the control to be optimal. 

As a special case of the problem above consider the situation where: 

• we have an infinite time interval, i.e. t 0 = O.tj —* oc 

• we have quasilinear dynamics, i.e. equation (3.3) can be put in the form 

x = A(x)x + B(x)u (3.10) 


and the cost function is quadratic, and does not explicitly depend on time, i.e.: 


£(x. u. r) = £(x, ii) = ^ (x T Qx + u T Ruj 


(3.11) 


Under these conditions we get the following expression for the optimal control: 


u op t = -R B (x) 


loTw / dJ (xY 


dx 


(3.12) 


while H-J-B-equation simplifies to: 


dx 


o = x t Qx + 2 ^ Ai x)x- B(xf dJ(x) T 


dx 


dx 


(3.13) 


Remarks: 


• We have set — 0 in equation (3.7) because we have an infinite time 

problem. 

• We see that all we need to construct the optimal control is to find an expression 
for ; we do not explicitly require an expression for J(x). 
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• It can be shown that given stabilizability and detectability assumptions the 
closed loop system will be stable [45]. 


Now assume that the infinite time, nonlinear, quadratic cost optimization problem is 
sufficiently well behaved for to be continuously differentiable (this assumption 

is often made in deriving the H-J-B equation). Then according to Section (2.1) there 
exists a matrix P(x) such that: 


a/(x)\ T 

dx ) 


P(x)x 


(3.14) 


Using this representation for 


we are led to consider the equation: 


0 = x T Qx + 2x T P T (x)A{x)x-x T P T {x)B{x)R- 1 B T (x)P(x)x (3.15) 

= x T (Q + P r (x)y4(x) + /l T (x)P(x) — P t (x)B(x)R~ 1 B t (x)P(x )) x(3.16) 


One solution to (3.16) can be found by setting P(x) = P(x) where P(x) is the positive 
definite symmetric solution to the matrix Ricatti Equation: 


0 = Q + P{x)A(x) + A T (x)P(x) - P(x)B(x)R 1 B T (x)P(x) (3.17) 


Let us examine the solution we have obtained in more detail. Although: 

P{x)x (3.18) 

satisfies equation (3.16), it has to meet certain requirements before we can say that 
we have a solution to the H-J-B-equation (3.13) (see also [45] for conditions that 
guarantee that a control is optimal). The requirements are: 

(1) P(x)x must be the gradient of some scalar function. That is, there must be 
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some scalar function V’(x) such that: 


(^) T = F(x)x (3.19) 

At this point we have no guarantee that there exists such a scalar function. 

( 2 ) If it is true that P(x)x is the gradient of some scalar function, then V'(x) has 
to furthermore satisfy the conditions (3. 8 ), (3. 9). 

The following lemma shows that if that P is the positive definite solution to (3.17) 
the second condition will be met. 


Lemma 3.2.1 If 



P(x)x 


(3.20) 


where P(x) is a positive definite matrix for all x and V r ( 0 ) = 0 then we will have: 


V(x) >0 Vx 7 ^ 0 


(3.21) 


Proof: 

Fix x and consider V(px) as a function of the scalar variable p. Then 


/o' 

(3.22) 

p 

(dV\ 


L 

(ajJ ^ 

\ / fLX 

(3.23) 

/:< 

px T P(px)j xdp 

(3.24) 

f /iX X Amin [^(/^x)] d/J, 

Jo 

(3.25) 


where A m i n [P(px)J denotes the minimum eigenvalue of P(px). Since P(x) is positive 
definite for all x we have A^n [P(/rx)j > 0 , and it follows that l/(x) > 0 for all x^O. 

Q.E.D. 
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With Lemma 3.2.1 in hand we can thus conclude that 


P(x)x (3.26) 

with P(x) the positive definite solution to equation (3.17), will be a solution to the 
H-J-B equation for those cases where P(x)x actually is the gradient of some scalar 
function. (Note that when we have a linear system, P does not depend on the state 
and thus P(x)x = Px will satisfy the gradient condition). 


3.2.2 The Continuous Ricatti Design Feedback Law 

For the continuous Ricatti design method we assume that the nonlinear system has 
been put in the quasilinear form of equation (3.10). We then use the form of the 
solution shown in the previous section without regard to whether condition (3.19) is 
satisfied or not. Specifically we use the following state feedback law: 

u = — A'(x)x (3.27) 

where 

A'(x) = R- l B T {x)P(x) (3.28) 

and P(x) is the symmetric positive semi-definite matrix solution to the matrix Ricatti 
equation: 


0 = Q + P(x)A(x) + A t {x)P(x) - P(x)B(x)R~ 1 B t (x)P{x) (3.29) 


Remarks: 

• The method we propose is generic in nature since a large class of nonlinear 
systems can be put in quasilinear form (see Section 2.1). 
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• We have found that when using the feedback law (3.27) the responses of the 
closed loop system could be influenced in a predictable manner by our choice 
of Q and R. 

• The feedback law is based on the dynamics of the system “frozen” at the current 
state. It “ignores” future behavior of A(x) and B(x), hence we refer to it as 
a control with zero look-ahead (see Chapter 4 for more discussion on the look- 
ahead concept). 

• Because condition (3.19) may not be satisfied we cannot claim optimality nor 
global stability for the closed loop system. However we have found that contin- 
uous Ricatti design controllers generally are better behaved than designs based 
on the linearized dynamics of the systems (see also Section 3.3.3 for further 
discussion of the stability properties). 


Example 3.2.1 Consider the nonlinear system: 


x = f(x) + g(x)tt (3.30) 


0.5 sin(x 2 )xi + x 2 


1 


(3.31 

+ 


u 

X\ — X\ X 2 


1 




One approach to developing a controller for this system would be to find a feedback 
law based on the linearized dynamics of the system. For example, linearizing ( 3.31) 
at the equilibrium point x = 0, gives: 


Sx 


df 

dx 


6x + 

x=0 


<X g«) 

du 


u 


x=0 


0 1 


Cr, ' 


1 




+ 


1 0 


— 1 

s 

1 


1 


(3.3-2) 

(3.33) 


A Linear Quadratic Regulator feedback law based on the linearized dynamics and 
cost function: 


J = 



(bx T Q6x -f uRu'j dt 


(3.34) 


53 



where: 



\l ol 



o 

II 

So 

II 

(3.35) 


L 0 ‘J 

gives the control law: 

U = —Kti n X 

(3.36) 

where: 

Kun = 1.366 1.366 

(3.37) 


If we apply this constant gain feedback to the full nonlinear system we get closed 
loop responses as shown in the phase plane plot of Figure 3-1. The closed loop 
system is locally stable (as it should be because its linearization is stable — see also 
Section 3.3.1) — but eventually becomes unstable. 

However if we apply the continuous Ricatti design method to the same system using 
the same Q and R, and the following quasilinear form for the dynamics of equa- 
tion (3.31): 

i\ 0.5sin(;T2) 1 .iq 1 

+ u (3.38) 

X 2 1 — X 2 1 

we get the “non-local” stable behavior shown in the phase plane plot of Figure 3-2. 
We also see that for the continuous Ricatti design method, the trajectories starting 
in the quadrant where X\ < 0,x 2 < 0 take a more direct route to the origin. 

3.2.3 Controllers for Combined Tracking and Regulation 

In this section we develop methods to achieve our second goal for the nonlinear design 
methodology, i.e. to find design variables that we can use in the synthesis process 
to adjust system responses. Our approach is to construct an appropriate quadratic 
regulator problem and then use the feedback law (3.27) given by the candidate so- 
lution to the H-J-B equation discussed in the previous section. We have found that 
despite the fact that the candidate solution may be non-optimal/unstable we could 


54 




Figure 3-2: Phase Plane Plot for Continuous Ricatti Design 
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still modify system behavior in a predictable way by adjusting the cost function (the 
Q and R matrices in equation (3.11)). 

Our goal is to find a control law that causes output variables, say y = Cx, to 
track desired values, say y^, and at the same time cause other state variables, say 
77 , to remain small. For example, we can think of y as variables that describe the 
gross motion of a flexible manipulator and 77 as its flex mode deflections. Of course 
depending on the system we are dealing with, the physics of the problem may not 
allow us to simultaneously achieve both accurate tracking (y % y^) as well as good 
regulation of the 77 variables. However by using a cost function approach, we have 
been able to conveniently make the trade-off between tracking accuracy and amount 
of regulation. 

The generality of this combined tracking and regulating goal becomes clearer if we 
recall (see [59] and Section 1 . 1 . 2 ) that through partial feedback linearization we can 
generally transform multivariable non-linear systems of the form: 


to the following form: 


x = f(x) + B(x)u 
= A(x)x + B(x) u 
y = Cx 


' ^ ll ' 


a x (x) 


' /3f(x)u ' 

y { 2 ri) 

— 

a 2 (x) 

+ 

0 a(x)u 

7 #( r m) 

ovn 


&m (x) 


. /£(x)u 


def 


Qj(x) 

«2(X) 


+ E(x) U 


ttm(x) 


(3.39) 

(3.40) 

(3.41) 


(3.42) 


(3.43) 
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V 


w(*7, x, u) 


(3.44) 


where y\ r) denotes the r’th derivative of y t . Here we have a generic separation of our 
state variables into two classes viz. 

• y = [yif ■ Vm ] T which are output variables which must track the desired tra- 
jectories y d(t) 

• r) which are other state variables w T hich we may want to keep small. 

Note that the continuous Ricatti design method does not require us to transform our 
dynamics to the partially feedback linearized form shown here — the analysis above 
only serves to show that the optimization problem is generic in some sense. 

Let us now construct the optimization problem for our quasilinear system. We will 
structure the problem in such a way that we remain with a regulator type problem 
(as opposed to an explicit tracking type problem [28]). Therefore we assume that the 
variables y d are the outputs of some reference model, that will generally be a stable 
linear system driven by an external input, viz: 

k d = A d x d + B d r (3.45) 

Yd = C d x d (3.46) 


where: 

x d is the / dimensional state vector of the reference model, 
r is an external input driving the reference model. 

We can now augment our basic plant dynamics (3.40) with our reference model to 
form the system: 


z 


F( z)z + G(z)u + Tr 

C -C d 


e 
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z 


(3.47) 

(3.48) 



where: 


x 


def 

z = 


Xd 


def 

e = y - y d 


(3.49) 


and 


F{ ■) = 


A(x) 0 

0 A, 


G(z) = 


5(x) 

0 


r = 


o 

5 , 


(3.50) 


In order to define a regulator problem we will temporarily set r = 0 but will show in 
Lemma 3.2.2 that under certain conditions we can get y to track y d even when r ^ 0. 


The cost function for our problem will penalize: the tracking error e = y — y d, 
excursions of those states which we want to regulate, as well as the control input, viz: 


J = J (e T Q y e + p x x T Q x x + p 2 n T uj dt (3.51) 

= j ((y - yd) T Q y {y - yd) + fj, x x. T Q x x. + p 2 u r u) dt (3.52) 

= f J (z t Qz + u T Ruj dt (3.53) 


where: 



C T Q y C + p x Q x -C T Q y C d 
-CjQyC CjQyCd 


(3.54) 


With the cost function as defined above (and r = 0) we have a problem in a format to 
which we can apply the continuous Ricatti design method as discussed in Section 3.2.2. 
We will furthermore show that independent of the question whether the continuous 
Ricatti design method is optimal or not we can get y to track y^ under certain 
conditions. To do this we need the following result by Kwakernaak and Sivan [33]: 


Theorem 3.2.1 Consider the time-invariant stabilizable and detectable linear sys- 
tem: 


i = FfZ + Gju 

y = c f z 
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(3.55) 

(3.56) 


where G / and C / are assumed to have full rank. Consider also the cost function: 

fOO . 

j ( z t CJQCjz + u r i?uj dt (3.57) 

with Q , i? positive-definite and symmetric . Let 

R = P 2 N (3.58) 

with N positive definite and p a positive scalar. Let P p be the steady-state solution of 
the Ricatti- equation: 


-PM) = CjQC, + PM)F / + FfPM)-PM)G,R-'G^PM) (3.59) 

PMC = 0 (3.60) 

Then the following facts hold: 

1. The limit: 

lim P p = P 0 (3.61) 

exists . 

2. Let z p (t),t > 0, denote the response of the controlled variable for the regulator 
that is steady-state optimal for R — p 2 N . Then 

TOO 

lrmy^ dt = z T (0)P o z(0) (3.62) 

(i.e. the term / 0 °° u r A'udt remains finite as p — ► OJ 

3. If dim(y) > dim(u) then Po ^ 0 

4 ■ //dim(y) = dim(u) and the numerator polynomial 't(.s) of the open-loop trans- 
fer matrix H(s) = Cf(sl - F f )~ l G f is nonzero, P 0 = 0 if and only if all the. 
zeros of the numerator polynomial [33] have non-positive real parts. 
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5. If dim(y) < dim(u), then P 0 = 0 if there exists a rectangular matrix M such that 
the numerator polynomial 'P of the square transfer matrix Cj(sl — A)~ l GjM is 
non-zero and has zeros with non-positive real parts only. 

6. If dim(y) = dim(u) and P 0 = 0 then: 

N 2 (hmG T f ^j = WQiCf (3.63) 

where W is some unitary matrix. 

We can now show when we can get y — > ► independent of the question of optimality 

(see also [53] for related results in a Loop Transfer Recovery setting). 

Theorem 3*2.2 For the system (3,47) } (3.48) and cost-function (3,53) assume that: 

• G(z), [C — Cd ] are full rank for all z with rank(G) = rank ([C — Cd]) 

• Q y = I ,R = p 2 1 and p is positive 

• (F(z),G(z)} is stabilizable for all z 

• {[C — Cot], F( z)} is detectable for all z 

• The transfer function [C — Cd](sl — F(z)) _1 G(z) is minimum-phase for all z 
If we apply the control: 

u = -A"(z)z (3.64) 

where: 

I<{ z) = \g{z) t P{z) (3.65) 

P 2 

and P(z) is the unique positive semi-definite solution to: 

0 = Q + P(z)F(z) + F t (z)P(z) - fp(z)G(z)G T (z)P(z) (3.66) 

P 2 
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with Q as in equation (3.54) then if (i x — 0 we will get 


y yd 


(3.67; 


as p — > 0. 

Proof: 

If we apply the control (3.64) to the system (3,47) we get for the closed loop dynamics: 


z = 


F(z) - -irG(z)G r (z)P(z) 


z + Tr 


1 


pz = pF(z)z — G(z) (z)P(z)J z + />rr 


Now taking the limit as p — ► 0 we get: 


1 


0 < G(z) ^-G J (z)P(z)J z 


By using Theorem (3.2.1) we get that as p — ► 0: 


(3.68) 

(3.69) 


(3.70) 


0 <- — G(z)B r (z) 


C-C d 


(3.71) 


where W( z) is a unitary matrix. Premultiplying by G' r (z) we get: 


0 G T (z)G(z)W(z) 


C - Cd 


(3.72) 


and using the facts that G T G is non-singular because G is full rank, and W is unitary 
we get: 

0 <— Gz - C d z - y - y d (3.73) 


.Q.E.D. 


Remarks: 


61 



• The above theorem shows that we asymptotically get perfect tracking when 
fj, x = 0 and p — > 0. By increasing p x we have been able to trade off tracking 
accuracy relative to the amount of regulation on the state variables which we 
want to keep small. 

• Although we may get perfect tracking with p x = 0 we are still not guaranteed 
stability as can easily be seen by examining equations (3.43) and (3.44). For 
instance although y may be well behaved rj may be unstable. 

• To obtain perfect tracking through partial feedback linearization one typically 
requires that E(x) in equation (3.43) be invertible. However we see that by using 
the continuous Ricatti Design method we do not need to perform the partial 
feedback linearization step and thus do not require that E(x) be invertible. 

• We have seen that we asymptotically achieve perfect tracking as p — > 0. This 
may give rise to the concern that we will get large control inputs. Fortunately 
we note that the control input requirements are essentially determined by the 
response/bandwidth of the reference model. Hence we can ensure that we have 
reasonable control inputs by suitably limiting the bandwidth of our reference 
model. 

• If we want proportional plus integral control we can further augment our system 
with states which are the integral of the tracking error. This will result in the 
augmented system 

z = F(z)z + (j(z)u + Tr (3.74) 

where now: 

z 

e/ 

F( z) 
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= [x, x rf , e/] T 

= y-y d 

A(x) 0 0 

= 0 A d 0 G{ z) 

C -C d 0 


(3.75) 

(3.76) 

B(x) 0 

0 r = B d (3.77) 
0 0 



To see how this results in a system with P-I-control let us examine the resulting 
control law. The closed loop control will be of the form: 


u = —K( z)z 


def 



K d 


Kj 


x 

Xd 

e/ 


(3.78) 


(3.79) 


(where we have suppressed the dependence on z of K x , K d , A'/ for notational 
convenience). If we assume that [C — C d \ is full rank we can find an invertible 
matrix T : 


ji def 


c -Cd 

ll 1 T22 

which we can use to transform to a new set of state variables, viz: 


(3.80) 


e 


= T 


x 


x r 


Xd 


(3.81) 


We can now rewrite the control law in a form that exhibits the proportional 
and integral feedback terms, as follows: 


u = — 

A' x K d 

r _1 r 

X 


^ ■* 


Xd 


Kiei 


def 



- A'/e/ 


Ap rop e h, nt ei A r x r 


(3.82) 


(3.83) 

(3.84) 


Figure 3-3 gives a block-diagram of this control law. (Note that for the quasi- 
linear system K pTop , K, nt , K r all depend on z although we have not explicitly 
shown this dependence). 

• The approach we used in this section also provides a useful solution to the 
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BLOCK-DIAGRAM AFTER 

aJ*H transformation 



Figure 3-3: Block Diagram of Augmented System with Integral Control 
so-called LQ-Servo problem [4] for linear time invariant systems. 


3.3 Stability Analysis 

In this section we examine the stability properties of the continuous Ricatti design 
method. It has been our experience that the method results in closed loop systems 
with good stability properties if we appropriately adjust the cost function. Our 
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approach will be to motivate our expectation of stability through analytical means 
but accept that we will have to assess the actual stability of the closed loop systems 
through computational methods. The computational methods will be discussed in 
Section 3.4, while we consider the following issues in this Section: 

• Local stability 

• Existence of multiple equilibrium points 

• “Non-local"’ stability 

We will prove that the control law will result in closed loop systems which are locally 
stable in the sense of Lyapunov (i.s.L.) and that the closed loop system will generally 
have a single equilibrium point at the origin. By using results from linear time varying 
systems theory we will furthermore show that the stability region can generally be 
expected to be large, but will not prove global asymptotic stability. 

We will use the following stability definitions in our analysis [59], [63]: 

Definition 3.3.1 Assume that the system: 

x = f(x) (3.85) 

has an equilibrium point at x = 0, i.e. f(0) = 0. Then 

• The equilibrium 0 is stable i.s.L. if for each e > 0 there exists a 8 — 6(e) > 0 
such that 

||x(0)|| < 6(e) => ||x(t) || < e, Vt > 0 (3.86) 

Otherwise the equilibrium point 0 is unstable 

• The equilibrium 0 is asymptotically stable i.s.L. if it is stable i.s.L. and if 
in addition there exists some q > 0 such that 

||x(0)|| < q => x(t) — * 0 as t — ► oc (3.87) 
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• The equilibrium 0 is exponentially stable if there exist constants r, a, A > 0 
such that 

||x(t)|| < a||x(0)||e _At , Vf > 0, V||x(0)|| < r (3.88) 

• If the equilibrium 0 is stable i.s.L., or asymptotically stable, or exponentially 
stable for all initial states, the equilibrium point is said to be globally stable i.s.L,, 
or globally asymptotically stable, or globally exponentially stable respectively. 

Note that exponential stability is the strongest form of stability since it implies asymp- 
totic stability as well as stability i.s.L. 


3.3.1 Local Stability Properties 

As a first step in analyzing the stability properties of the continuous Ricatti design 
method we prove that the method yields closed loop systems that are locally stable. 

To show local stability for the continuous Ricatti design method we use Lyapunov’s 
linearization method which is based on the following theorem [59]: 


Theorem 3.3.1 Consider the system: 

(3.89) 

(3.90) 


x = f(x) 
f(0) = 0 


Suppose that f(x) is continuously differentiable and let J( 0) be the Jacobian of f(x) 
evaluated at x = 0, i.e.: 


J{ 0) d = 


df 

<9x 


x=0 


(3.91) 


Under these conditions: 


• If all eigenvalues of J( 0) are strictly in the left half complex plane then 0 is an 
exponentially stable equilibrium point of (3.90). 
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• If one or more eigenvalues of J( 0) is strictly in the right-half complex plane 
then 0 is an unstable equilibrium point. 

• If all the eigenvalues of J(0) are in the left-hand complex plane with at least 
one eigenvalue on the juf-axis then one cannot conclude stability nor instability 
of the nonlinear system (3.90) from its linearized approximation. 

Theorem (3.3.1) makes it clear that in order to establish local stability of the contin- 
uous Ricatti design method it is sufficient to examine the linearized approximation of 
the closed loop system. With this in mind we next show how we can easily obtain the 
linearized approximation of a nonlinear system if we have its dynamics in quasilinear 
form. (Note that the proof below holds for any quasilinear realization of f(x), not 
only those obtained via equation (2.6)). 

Lemma 3.3*1 Consider the nonlinear system 

x = f (x) (3.92) 

= A(x)x (3.93) 

with f(x) continuously differentiable in an open region around 0. Then the Jacobian 
off(x ) is given by: 

= A{ 0) (3.94) 

x=0 


af 

dx 


Proof: 

Recall (see [46]) that the Jacobian of f(x) evaluated at Xo is equal to the matrix J(x 0 ) 
such that: 

f(XO + h) -|f^ ) - J(X ° )h -0 as h — * 0 (3.9.5) 

ll h ll 

It can furthermore be shown that the matrix J(x 0 ) is unique [46]. For our case we 
want the Jacobian of: 


f(x) = .4(x)x 


(3.96) 
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evaluated at x = 0. Now note that: 


f(h) -f(0) - 4(0)h 

INI 


||f(h)-f(0)-4(0)hj| 

INI 


4(h)h -0 - 4(0)h 

INI 

(4(h)- 4(0)) h 

INI 

< ||A(h)- 4(0)11, -||h|| 

INI 

= ||A(h)- A(0)|| t Vh^O 


(3.97) 

(3.98) 

(3.99) 
(3.100) 


so that: 


f(h) — f(0) — 4(0)h 

INI 


as 


(3.101) 


Q.E.D. 


We are now in a position to prove our assertion that the continuous Ricatti design 
method will yield locally stable controllers, viz: 

Lemma 3.3.2 Given the system: 

x = 4(x)x + B(x)u (3.102) 

and feedback law 

u = -K(x)x (3.103) 

where 

K(x) = R~ 1 B t (x)P(x)x (3.104) 

and P(x) is the positive-semi- definite solution to: 

o = Q + P(x)4(x) + 4(x) r P(x) - P(x)B(x)R~ 1 B t (x)P{x) (3.105) 

Assume 4(x),J5(x) is stabilizable and 4(x),Qz is detectable. Then the closed loop 
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system will be locally exponentially stable. 


Proof: 

The closed loop dynamics of the system is given by: 

x = f e /(x) (3.106) 

= A c ,(x)x (3.107) 

where: 

Ad(x) = .4(x) - B(x)K(x) (3.108) 

Since we have .4(0), 5(0) stabilizable and A(0),Q* detectable we know from the 
properties of Linear Quadratic Regulators that A c /(0) will have all its poles strictly 
in the left half plane (see e.g. [32]). Furthermore from Lemma (3.3.1) we know that 
the Jacobian of f c ;(x) is given by A c /(0). Local exponential stability of the continuous 
Ricatti design method then follows from Theorem 3.3.1. 

Q.E.D. 

At this point we have rigorously established the local stability of the continuous 
Ricatti design method. In Section 3.3.3 we will further analyze the stability of the 
method and show why we expect the stability region for the continuous Ricatti design 
method to be large. 

3.3.2 Single Equilibrium Point at the Origin 

If we applied an arbitrary control law to the system (3.102) we might get multiple 
equilibrium points for the closed loop system which would in general be undesirable. 
However the following lemma shows that if we apply the continuous Ricatti design 
method we will have only a single equilibrium point at the origin. 
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Lemma 3.3.3 Given the system, feedback law and stabilizability and detectability 
assumptions of Lemma 3.3.2, the closed loop system will have a single equilibrium 
point at x = 0 

Proof: 

The closed loop dynamics are given by equations (3. 107), (3. 108). Clearly at x = 0 we 
will have x = 0. It remains to be shown that x/0, Vx 7^ 0. Since j.4(x), B(x),Q*} 
is stabilizable and detectable and K (x) is the feedback gain from the steady state 
regulator problem, the matrix A c /(x) = (A(x) — Bfx)Kfx)) will have all its eigenval- 
ues in the open left half plane. Thus .4 c ;(x) will have an empty nullspace and thus 
A c i(x)x / 0 Vx ^ 0. 

Q.E.D. 


3.3.3 Non-local Stability Analysis via Linear Time Varying 
Regulator Theory 

The continuous Ricatti design method outlined in Section 3.1 can also be motivated 
from a stability standpoint by using results that are available for linear time varying 
systems. Our goal will be to provide a rationale of why we expect the continuous 
Ricatti design method to yield closed loop systems with large stability domains. The 
connection between (3.102) and linear time varying systems comes by noting that 
along a given trajectory x(t) the matrices A(x), B(x) become functions of time. So 
we can think of our quasilinear system as an ensemble of time varying systems "p a ~ 
rameterized” by the trajectory x(t), viz: 

x = A(x(t))x + B(x(t)) u (3.109) 

« A(t)x + B(t) u (3.110) 
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Chapter 4 will discuss a method which further exploits the analogy with linear time- 
varying systems. We should point out that because we are not dealing with a true 
linear time varying system the analysis we provide here does not guarantee stability, 
but only gives us some expectation of stable behavior. 

In order to motivate the stability of the continuous Ricatti design method we need 
two results for linear time varying systems (LTV systems). The first result we need 
comes from Kokotovic et al [30]. Consider the following optimization problem for 
LTV systems: 


For the system: 


x = A(t)x + B(t ) u with x(0) = x 0 


find the control which minimizes: 


J = ~ j [x T C T Cx + u T Ruj dt 


(3.111) 


(3.112) 


subject to: 


x{T) = XT 


(3.113) 


Reference [30] shows that if T is large with respect to the system dynamics and 
certain conditions outlined below are met, then the solution to the problem above 
exhibits "boundary layers”. That is, the full optimal control can be approximated by 
the solution of two simpler problems, one at t = 0 and one at t = T. The following 
theorem slightly adapted from [30] makes these ideas more precise: 

Theorem 3.3.2 For the system (3.111) and cost-function (3.112), assume that: 

• A(t ), B(t) are continuous for all t £ [0, T] 
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• The eigenvalues of the Hamiltonian matrix: 


A{t) -B{t)R~'B{t) 


-C T C 


-A r (t) 


(3.114) 


all lie off the imaginary axis. 

• The pair A(t), B(t) satisfies the frozen controllability condition, that is: 


rank B(t), A(t)B(t), . . . , A n 1 (t)B(t) = n \/t € [0, T] 


(3.115) 


where n = dim(x) of the system. 

Then there exists an t\ > 0 such that for all e = f ^ £ (0, ej] the optimal state trajectory 
and optimal control satisfies: 


x(f) = xfit) +x r (<r) + 0(e) 
u (t) = ufit) + u r (<) + 0(e) 


(3.116) 

(3.117) 


where the initial and terminal boundary layer controls are given by: 


u fit) = ~R- 1 B t ( 0)A,(<) 
u r(t) = -R- 1 B T {T)\ t (<t) 


(3.118) 

(3.119) 


def rp , 

<7 — T — t 


A fit) = P(0)xi(t) 

\ r (a) t: f N(T)x r (a) 


(3.120) 

(3.121) 

(3.122) 
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and 


^ = [A(0) - B{0)R 1 5 r (0)P(0)] x, x,| <=0 = x 0 (3.123) 

^ = -{a(T)-B(T)R-'B t (T)N(T)\x t x r |„ 0 = x 7 (3.124) 

and .P(O) is the positive definite solution of the algebraic Ricatti equation: 

0 = C T C + V(t)A(t) + A(t) T V(t) - r(t)B(t)R- l B T (t)V(t) (3.125) 

at t = 0, and N(T) is the negative definite solution to the same equation at t — T . 

The usefulness of this theorem is that it enables one to find an asymptotic approxima- 
tion to the solution of the time-varying optimal control problem by solving two time 
invariant Linear Quadratic Regulator problems, one at t = 0 and one at / = T. The 
boundary layer notion can be explained as follows. Under the frozen controllability 
assumption (3.115) the dynamics of both equation (3.123) and equation (3.124) will 
be stable [40]. Hence for T 1 we will have: 

x,| t=r « 0 (3.126) 

Xr| a=T ~ 0 (3.127) 

Thus at t « 0: 

X ~ X; (3.128) 

u % -R- 1 5 r (0)P(0)x ; (3.129) 

That is, the initial (boundary layer) behavior of the closed loop system is dominated 

by the dynamics of equation (3.123). Similarly at t ss T => a ss 0 

x « x r (3.130) 

u « -R- 1 B T {T)N(T)x r (3.131) 
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So we see that the terminal (boundary layer) behavior of the closed loop system 
is dominated by the dynamics of equation (3.124). From equation (3.129) it also 
becomes clear that to approximately compute the initial control for the given problem 
we do not need to know the future values of A(£), B(t) provided T — > oc. (Note: If 
Xj = 0 then x r = 0 Vt and the approximation will be even better). Furthermore 
the result implies that the control input (3.129) will be asymptotically correct for an 
ensemble of linear time varying systems that have the same value for .4(0). B( 0). 

The second result we need comes from receding horizon control theory (see Chapter 4 
for more detailed discussion and references). Receding Horizon Control uses the 
following strategy: 

1. At each time instant t solve a linear quadratic regulator problem over a horizon 
starting at the current time t and ending at the time t + T subject to the 
terminal constraint x(t + T) = 0. 

2. Then at each time instant t apply the initial control, i.e. u (t), found from the 
regulator problem that we have just solved. 

The following theorem slightly adapted from [36] summarizes the fundamental sta- 
bility result for receding horizon control. 

Theorem 3.3.3 Let u‘(r),r 6 [ t,t-\-T ] be the control time history which minimizes 
the cost: 

J = ^ J (x. T C T Cx + u T Ruj dr (3.132) 

for the system (3.111) subject to the terminal constraint x(t + T) = 0. Assume that: 

• R is positive definite 

• Q C T C is positive semi-definite 

• the pair A(t),B(t) is uniformly completely controllable for some T c > 0 
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Then 


1. If we choose T > T c and at each time t we apply the control u = u *(t) then the 
closed loop system will be uniformly asymptotically stable. 

2. The control u *(t) is given by: 

u(t) = -R- 1 B T (t)P(tJ + T)x(t) (3.133) 

where P(t,t + T) = M~ l (t,t T T ) and M{t,t -f T) is found by integrating the 
time-varying Ricatti Equation: 

8M(t , t + T) . , „ 

^ " = B(t)R~ x B { t)- M{tA + T)A t (t)~ A{T)M{ T ,t + T) 

—M(r,t A T)C T CM(r,t + T) (3.134) 

backwards from r = t + T to r = t subject to: 

M(t + Tj + T) = 0 (3.135) 

Note that the receding horizon control method is no longer optimal for the original 
cost function (3.132) (see [35] for a bound on the cost incurred by the receding horizon 
control law). 

By combining Theorems (3.3.2) and (3.3.3) we can now provide a stability argument 
for the case where we are applying the continuous Ricatti design method to linear 
time-varying systems, viz: 

• From Theorem (3.3.3) we know that the system (3.111) will be stable if we 
apply the receding horizon control described above. 

• In stead of finding the control (3.133) by integrating the matrix Ricatti differ- 
ential equation (3.134) over the interval [t. t + T] we can find an asymptotically 
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correct value (as T — * oo) for u’(£) by using Theorem (3.3.2). The approximate 
value of u*(<) is given by: 


u{t) = -R~ 1 B T (t)P{t)x{t) (3.136) 

where P(t) is the positive definite solution to: 

0 = C T C + P(t)A(t) + A T (t)P(t) - P{t)B(t)R~ 1 B T (t)P(t) (3.137) 

So we expect stable behavior for linear time-varying systems insofar as the 
asymptotic approximation is valid. 

We see that the procedure outlined above amounts to the continuous Ricatti design 
method as applied to Linear Time Varying Systems. To motivate our expectation 
of closed loop stability for the case where we apply the continuous Ricatti design 
method to nonlinear systems, we treat the quasilinear system (3.109) as if it were a 
linear time varying system along a given solution trajectory as in equation (3.110). 
Here we exploit the fact that to find the asymptotically correct value of the receding 
horizon control input (u* in equation (3.133)) we only need to know the current values 
of A(x(J)), B(x(t)), However note that while we can generally assume that A(t ), B(t) 
are well-behaved for the true linear time-varying case, we have no such guarantees 
for A(x(t)), B(x(t)) and hence cannot guarantee stability for the case when we are 
dealing with nonlinear systems. 


3.4 Stability Region Assessment 

In the previous sections we have analyzed the stability properties of the continuous 
Ricatti design method and have shown rigorously that the closed loop systems will be 
locally stable. We also indicated why we expect the stability domains to be relatively 
large. In this section we will examine computational methods to assess the size of 
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the stability domain. The first method will provide a straightforward procedure to 
assess a lower bound on the size of the stability domain while the second method will 
provide a computationally more complex procedure to determine the exact extent of 
the stability domain. Although not discussed here further we should point out that 
simulating the closed loop dynamics still provides a very useful method for assessing 
the stability of nonlinear systems. 


3.4.1 Stability Domain Assessment Using Quadratic Forms 

The first method we use to assess the stability domain of the nonlinear closed loop 
system will take a Lyapunov function approach and will provide a lower bound on 
the size of the stability domain (see also [21], [64] for related results). In order to 
develop the method we will use La Salle’s theorem [37] which is stated below. Before 
stating La Salle’s theorem we will need the following definitions (see [59]): 

Definition 3.4.1 1. A scalar continuous function V(x) is said to be locally pos- 

itive definite ifV(O) = 0 and in an open neighborhood of the origin: 

x / 0 => !/(x) > 0 (3.138) 

2. IfV( 0) = 0 and the above property holds over the whole state space, then V'(x) 
is said to be globally positive definite. 

3. A scalar continuous function V(x) is said to be negative definite if—V(x) is 
positive definite. 

Theorem 3.4.1 (La SalleJ Consider an autonomous system of the form: 

x = f(x) (3.139) 
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with f(x) continuous and f(0) = 0. Let V'(x) be a locally positive definite function 
with continuous first partial derivatives and suppose that: 

• for some l > 0 the set : 

ft/ = f {x)V’(x) < /} (3.140) 

is bounded . 

• V is bounded below on ft/ 

• V{x) <0 VxGfi, 

• the set: 

S = f {x € ftj|V(x) = 0} (3.141) 

contains no trajectories of (3.139) other than the trivial trajectory x(t) = 0. 

Then the equilibrium point 0 of (3.139) is asymptotically stable. 

The above theorem provides us with sufficient conditions that guarantee that trajec- 
tories starting in a domain ft/ will tend to 0 as t — ► oc provided that we can find a 
suitable function V^x). The following corollary shows how we can always construct a 
function V (x) that satisfies the conditions of La Salle’s theorem and provides us with 
a well defined region ft/ if the linearized version of (3.139) is locally asymptotically 
stable (note that Lyapunov’s linearization method, Theorem 3.3.1, only states that 
some stable neighborhood exists but does not show how to determine it). 

Corollary 3.4.1 Consider a system: 


X = f (x) 


(3.142) 


with f(x) continuously differentiable and f(0) = 0. Let 


A 


def 


df 


dx 


x=0 


(3.143) 
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Assume that A has all its eigenvalues in the open left half plane and let P 0 be the 
positive definite solution to the matrix Lyapunov equation: 


PqA + A t P q — —I (3.144) 

(Note: Pq > 0 is guaranteed to exist, since A is strictly stable see [59]). Define: 

V(x) = f ^x r F 0 x (3.145) 

S ~ = fa {(x) = ° } (3.146) 

and let l be the lower bound of the values of V / (x) on the set S with the exception 
of the origin (if S contains only the origin set l = oc). Then the origin will be 
asymptotically stable for every solution originating in the region f lu defined by: 

Liu = f {x|V"(x) < l - s} (3.147) 


where £ > 0. 

Proof: 

We will prove this corollary by showing that when / is finite we satisfy the requirements 
of Theorem (3.4.1) and that when / is infinite, (3.142) will be globally asymptotically 
stable by application of Lyapunov’s direct method [59]. 

Before proceeding we need the following from [63]: 

• Claim: V(x) is a (locally) negative definite function 
Proof of Claim: 

We can write: 

f(x) = Ax + r(x) (3.148) 
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with r(0) = 0 since Ax is the first term in the Taylor expansion of f(x). Now: 


V 

= x T P 0 f(x) 

(3.149) 


= x T P 0 (Ax -|- r(x)) 

(3.150) 


= |x r (P 0 A + A t P 0 ) x + x T P 0 r(x) 

(3.151) 


= - ^;X T x + x r p 0 r(x) 

(3.152) 

Now 

x T Po r(x) < ||x||*||r(x)|| 

(3.153) 

where a is the maximum eigenvalue of P 0 . Since r(x) is the remainder term of 

the Taylor expansion of f(x) we know that: 



lim "f)" =0 

INI-o ||x|| 

(3.154) 

and hence we can find 

an open neighborhood of the origin say J\f such that: 


||r(x)|| < -^||x|| Vx G Af 

(3.155) 

Hence 

x T P 0 r{x) < h|x||||x|| VxeAT 
4 

(3.156) 

and subsequently 

v < -bc T x + h|x || 2 VxeAT 

£ T 

(3.157) 


so that V is locally negative definite around the origin. 


Let us now consider the case where / is finite. In this case Q[ c will be a finite hyper- 
ellipsoid centered at the origin since Pq is positive definite and thus: 


V'(x) = x T P 0 x < l Vx € 


(3.158) 
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Furthermore: 


V’(x) ^ 0 Vx e Qu, x ^ 0 (3.159) 

(otherwise we could find x* such that V 7 (x*) < l and V"(x*) = 0 which would con- 
tradict our assumption that / is a lower bound). Since V is negative definite in a 
neighborhood of the origin and V is continuous we see from the intermediate value 
theorem [46] that V < 0 Vx € ft;. Asymptotic stability then follows from Theo- 
rem (3.4.1) since S = (x £ Q e |V(x) = 0} contains only the origin. 

The remaining case we have to consider is when / is infinite, i.e. V / 0 Vx ^ 0. 
Since V is at least locally negative definite and V is continuous it then follows from 
the intermediate value theorem that V < 0 Vx ^ 0 . Global asymptotic stability 
then follows from Lyapunov’s direct method by noting that V(x) = x T P 0 x is radially 
unbounded [59]. 

Q.E.D. 


With Corollary (3.142) in place we can now specify a computational procedure which 
will give a lower bound on the size of the stability region for a closed loop nonlinear 
system, i.e. a system of the form (3.142). To compute an estimate of the stability 
region we do the following 


1. Compute: 


A = 


df 


dx 


x=0 


2. Solve the matrix Lyapunov equation: 


PoA + A t P 0 = -I 


(3.160) 


(3.161) 


3. Find the minimum value of C(x) such that C(x) ^ 0 and V'(x) = 0. Some 
possibilities for doing this are: 
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(a) Compute V and V for several points along ‘‘rays” emanating from the 
origin to find / = f “smallest value of V where V = 0” . 

(b) Use a single variable rootfinding algorithm along “rays” emanating from 
the origin to find the first location where V = 0 along the ray. Among 
these locations find the value of x that minimizes V 

(c) Use a multivariable rootfinding algorithm to find locations where V = 0. 
Among these locations find the value of x that minimizes V . 

4. Once we have found the desired minimum value of V, say V m ini we know from 
Corollary (3.142) that all the trajectories starting within the hyper-ellipsoid 
defined by x T P qX = U min will tend to 0 as t — > oo 

Example 3.4.1 As an example of using this approach consider the following example 
system from Hahn [21]: 


x = — x + 2 x 2 y (3.162) 

y = -y (3.163) 


For this example it is possible to exactly determine the size of the stability domain 
(see [21] and Section 3.4.2). It can be shown that all trajectories starting in the region 
defined by 


y < 

1 

X 

x > 0 

(3.164) 

y > 

l 

X 

x < 0 

(3.165) 

y € ( oo, oo), 

for 

x = 0 

(3.166) 


will tend to the origin as t — * oo. 


Let us now use the quadratic form method outlined above. We get: 


-■I 


(3.167) 


x=0 
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Figure 3-4: Stability Domain Using Quadratic Form 


-1 0 
0 -1 

while P 0 from equation (3.144) evaluates to: 



(3.168) 


(3.169) 


Note that for this example the contours of constant V are circles. Searching for 
locations where V = 0 we get the results shown in Figure 3-4. The stability domain 
we find using this method is significantly smaller than the actual stability domain 
because it is constrained to be an ellipsoid. 
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3.4.2 Stability Domain Assessment using Zubov’s Method 


In this section we will discuss Zubov’s method (see e.g. [21], [64]) which can be used 
to determine the exact domain of stability for nonlinear systems of the form 


x = f(x) (3.170) 

Zubov's method has two important advantages for stability analysis: 

1. It provides a deterministic method to find Lyapunov functions for non-linear 
systems, as opposed to the usual procedure of “guessing" candidate Lyapunov 
functions. 

2. It provides a method to exactly determine the extent of the stability domain. 
This is different to standard Lyapunov theory which tends to provide results 
that either: 

• guarantee the existence of some stable ^-neighborhood of an equilibrium 
point but give no indication of the extent of the stability domain. 

• guarantee that systems are globally stable, (finding Lyapunov functions 
for this case is notoriously difficult). 

The advantages come with a price however. It requires one to solve a partial differ- 
ential equation (P.D.E.). Zubov’s main theorem makes this more precise [64], [21]: 

Theorem 3.4.2 LetU be an open domain of the state space and let U be its closure; 
U shall contain the origin. Necessary and sufficient conditions for U to be the exact 
domain of attraction of the equilibrium of 

x = f(x) (f (0) = 0) (3.171) 

is the existence of two functions v(x) and <^>(x) with the following properties: 
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1. v(x) is defined and continuous in U , y(x) is defined and continuous in the whole 
space 

2 . y?(x) is positive definite for all x 

3. v(x) is negative definite , and for xeW,x/0 ( the inequality 0 < t?(x) < 1 
holds. 

U Y € — U, then lim x _ y v(x) = 1. furthermore = 1 provided that 

this limit process can be carried out for x £ It . 

5 . 

dv 

— f(x) = ^(x)(i ! (x)- 1) (3.172) 

The following slightly modified example from Hahn [21] shows how the theorem can 
be applied by analytically finding a solution to (3.172). 


Example 3.4.2 Consider the system: 


x = —x + 2 x 2 y 

y = -y 


(3.173) 

(3.174) 


For this example, Zubov’s P.D.E. (i.e. equation 3.172) becomes: 


dv dv 

-(-x + 2 x 2 y) + -(-?/) = f(x.y){v - 1) 


(3.175) 


If we set ip(x,y) to the positive definite function: 


w*,*) =(^y + (£) 


(3.176) 


we can find the following solution to equation (3.175): 


v(x,y) = 1 -exp (-i (I)* - i (i)’ (j-2 




(3.177 
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We see that the curve xy — 1 defines the stability boundary. 


Numerical Solution of Zubov’s Equation using Finite Differences 

The example above demonstrates Zubov’s method, but it also shows that in general 
we will have difficulty in obtaining exact analytical solutions. We therefore consider 
a numerical solution procedure to the example above. Several methods for solving 
P.D.E.’s exist, such as finite difference, finite elements etc. (see e.g. [1]). Perhaps 
the simplest approach to solving (3.172) is to use a finite difference method which we 
shall use below. 

The usual approach to solving the P.D.E. (3.172) by finite difference methods will be 
to represent the unknown function, v(x) in our case, by its values at a discrete grid 
of reference points. For example, for a second order system we could use: 

Xi — Xo + ih (3.178) 

Vj = Vo+jh (3.179) 

as shown in Figure (3-5). Such an approach however rapidly becomes unmanageable 
if we attempt to solve the P.D.E. over the whole domain of interest at once. For 
example, say w r e want to solve Zubov’s P.D.E. for an 8’th order system over a grid 
of 10 points per dimension. This means that we need to solve for the values of v at 
10 8 points which translates into the numerical challenge of solving a set of 10 8 linear 
equations. Therefore we propose to solve the P.D.E. for smaller domains at a time. 
A straightforward way to accomplish this is to solve the P.D.E. along narrow solution 
strips as shown in Figure (3-6). To assess the size of the stability domain we “point” 
these solution strips into all the directions of interest. Note that this approach has 
the added advantage that it is amenable to parallelization. 

Let us now develop a finite difference solution for Zubov’s equation in more detail. 
We discuss the method as applied to a second order system noting that the same 
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Figure 3-6: Solution Strips for Zubov’s P.D.E. 


87 


procedure can be used for higher order systems. 


Figure (3-6) shows a general solution strip that is rotated by an angle 6 from the 
nominal x, y-axes. Let: 

• x r denote the co-ordinates of a point in rotated co-ordinates 

• x denote the co-ordinates of the same point in non-rotated co-ordinates 

Zubov’s equation in rotated co-ordinates then becomes: 

(3.180) 

To evaluate x r we note that the transformation from non-rotated to rotated co- 
ordinates is a simple linear transformation, viz: 


dv 

dx r 


x r — ip(v — 1) 


x r = Cx 


(3.181) 


where C is the direction cosine matrix relating the two sets of axes. For our 2-D 
example we will have: 


x r 

, C = 

cos(0) sin(0) 

, x = 

X 

y T 


— sin(0) cos(0) 


y 


So that the system dynamics in rotated co-ordinates becomes: 


x r = Cx 

(3.183) 

= Cf(C -1 x r ) 

(3.184) 

d = f r (x r ) 

(3.185) 


At this point we are in a position to set up the finite-difference approximations for 
Zubov’s equation. We develop the approximations in rotated co-ordinates because 
of the simplification we obtain when the solution strip becomes a rectangular grid 
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aligned with the x r axis (see Figure (3-6)). Let us denote the point: 


x r 


Xq + ih 

. yT . 


_ Vo + J h 


on the grid, by the pair Similarly let: 


(3.186) 


v(x\y r ) = V lo , f r (x r ) 


r 

K 

Sb 

k 

k- 

1 


l 

k. _T 

i 

. fl{x T ,y r ) _ 


i 

i 


A.x\y r ) = 


[3.187 


Using this notation we get the following finite difference approximations: 


• For a point, say (i,j), interior to the grid: 


dv 

dx r 

dv 


Vi+ij - Vi-u 

2 h 

Vjj+i ~ Vij-i 

2 h 


(3.188) 

(3.189) 


• For points on the edges we use one-sided finite-difference approximations as 
appropriate. For example for a point, say (i,j), on the left edge away from the 
corners: 


dv 

dx r 


dv 



K+l.j - Kj 

h 

Vm - Vij-1 

2h 


(3.190) 

(3.191) 


Using these approximations, Zubov’s equation becomes a set of linear equations in the 
variables V t j,i = 1 . . . M,j = 1 . . . N, one equation for each grid point. For example 
at a point (i,j) interior to the grid, we get: 


V i+lJ - Vi. 


id 


2 h 


fl, + 


Vij+i - Vj. 


2h 


l ) = 


ViAVij - 1 ) 


[3.1921 
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which becomes a linear equation in the unknowns VJ.j, U\j+i, U+ij, viz: 

- - ^Vy., - VJyKj + K. J+ , + - -?,j (3.193) 


Similar results hold for other points at the edges of the grid. Assembling all these 
equations we get a set of linear equations which we have to solve: 


Av = — b 


(3.194) 


In equation ( 3.194): 


v — Vi'i Vj.2 ••• Kjv ^2,1 V2„v ••• Vjv/,1 ... Vm,n (3.195) 

k = 'r’l.l ••• ^2,1 ••• V^ 2 ,.V ••• ••• Vm,n (3.196) 

and each row in A corresponds to an equation such as (3.193). In general the matrix A 
will be sparse as can be seen by examining the row corresponding to equation (3.193): 

.0 ... 0 -% 0 ... 0 -% %- 0 ...0 % 0 ... 0 

(3.197) 

Exploiting the fact that A is sparse will enable us to further reduce our memory 
requirements [9]. 

Example 3.4.3 Consider again example 3.4.2 for which we know the exact solution. 
Using the finite difference method outlined above to solve Zubov’s P.D.E. for this case 
we get the results shown in figures (3-7) and (3-8). The plots compare the values of 
the finite difference solution and the exact solution along a solution strip — values 
are plotted in rotated co-ordinates, i.e. along the x r -axis of Figure 3-6. We used a 
solution strip as in Figure 3-6 with the following parameters: 


6 = 30° 


h = 0.05 


(3.198) 

(3.199) 



Comparison of finite difference and exact solutions 


1 


0.9 h 


Salid line - numerical solution _* 


0.8 




Dashed line - exact solution 



Figure 3-7: Exact vs. Finite Difference Solution to Zubov’s P.D.E. 

(3.200) 

(3.201) 


M = 56 
N = 10 


We used <p as in example 3.4.2 with: 


a = 2 


(3.202) 


Numerical Conditioning of the Solution to Zubov’s Equation 

The approach outlined above to solve Zubov’s P.D.E. has yielded promising results 
but further work is required to analyze the numerical conditioning of solutions. For 
instance in example 3.4.2 the exact solution to Zubov’s P.D.E. is not well behaved. 
We see from equation (3.177) that when xy < 1 we have 0 < v(x. y) < 1. However 
for xy — 1 + e with 0 < £ <C 1 we see that v(x,y) will be large and negative. This 
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Figure 3-8: Error between Exact and Finite Difference Solutions to Zubov’s P.D.E, 


discontinuity in the (true) solution to the P.D.E. will generally cause difficulty for 
numerical solution methods. 

A further issue is that the choice of ^p(x,y) also affects the behavior of the solution 
as explained in the following. 


Zubov’s equation can be written as: 


|^f(x) = v?(x)(v(x) - 1) 


(3.203) 


Along solution trajectories of x = f(x) this equation can be written as: 


- = S,(x)(,-l) 


(3.204) 


or in reverse time r: 


dv 

dr 


-¥>(x)(t> - 1) 


(3.205) 
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So along a specific trajectory x(r) in reverse time we get: 

= -<?(x(r))(u - 1) (3.206) 

which has a solution of the form: 

v = 1 + cexp“^^ dr c = constant (3.207) 

The influence of p (which is a positive definite function) now becomes clearer. For 
example if p is relatively small then v will tend to remain small and grow quickly 
at the stability boundary where it must attain the value v = 1. On the other hand 
if p is relatively large then v will quickly approach its asymptotic value and it will 
numerically not be clear whether we are close to the stability boundary. 

Example 3.4*4 Figures 3-9 and 3-10 show the results of redoing Example (3.4.3) 
where the only change has been in the parameter a of the function p(x,y) of equa- 
tion (3.176). Using a = 0.1 (i.e. ip relatively large) we get the results in Figure 3-9 
which confirm that v quickly approaches its asymptotic value. Conversely when 
a = 10 (i.e. p relatively small), we see in Figure (3-10) that v remains relatively 
small. Note that in both cases the finite difference method produced accurate results. 


3.5 Computational Issues 

The continuous Ricatti design method requires in principle that we compute a feed- 
back gain: 

K(x) = R~ 1 B t (x)P(x) (3.208) 

by first solving a steady state Ricatti equation: 

0 = Q + P(x)A(x) + A(x) t P(x) - P(x)B(x)R- 1 B t (x)P(x) (3.209) 
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Comparison of finite difference and exact solutions 


Solid tin© — numerical solution 


Dashed line - exact solution 


alfa - 0.1 , theta - 30 deg 


i 



Figure 3-9: Exact vs. Finite Difference Solution to Zubov’s P.D.E. — large 


Comparison of finite difference and exact solutions 



Figure 3-10: Exact vs. Finite Difference Solution to Zubov’s P.D.E. — small 
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for each state x. To implement the control law we can: 


• compute the feedback gains online as the system evolves along a state trajectory. 

For this case we have found the following method to be effective: 

1. Assume we have an "old” solution to the Ricatti equation, say P old . As the 
closed loop system progresses along a trajectory keep using the feedback 
gain: 

A'(x) = -R~ x B T {x)P old (3.210) 

and “continuously” re-compute the residual A, where: 

A(x) = Q + P m A(x) + A T {x)P old - P old B(x)R~ 1 B T (x)P old (3.211) 


2. When 


l|A(x) 

\\Pm\ 




(3.212) 


for some predefined Sr > 0 recompute the solution to the Ricatti Equation 
and update P 0 i d with the newly computed value of P(x). 


• Precompute the feedback gains. Such an approach is examined in [62] where 
the feedback gains K (x) are computed and then approximated as a sum of some 
basis functions, viz: 

p 

#»j(x) = ]T a, jk i k (x) (3.213) 

k=i 

Using basis functions to approximate the feedback gains has the advantage that 
it reduces the storage requirements. 


Several methods exist to numerically solve the steady state Ricatti Equation. The 
methods can be classified as: 


• Non-Iterative methods, such as eigenvector [39] and Schur Vector Methods [38]. 

• Iterative methods, such as Kleinman's-method [29]. 
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We will summarize some of these methods below and introduce a new iterative method 
which we will refer to as Collar and Jahn's method. In each case we will assess the 
computational burden of the method by counting the number of “FLOPS” (floating 
point operations) required to solve the Ricatti equation. For comparison purposes we 
will count a scalar add and multiply,: 


a + (6*c) (3.214) 

as one FLOP, and count FLOPS for other numerical procedures according to the 
following table (see e.g. [50], [38], [5]). 


Operation 

Operation Count 
(FLOPS) 

Scalar add and multiply a + (b* c) 

1 

Matrix addition:/! + B A, B € R nxn 

n 2 

Vector dotproduct: x r y x, y € /? nxl 

n 

Matrix product: A * B A, B € R nxn 

n 3 

LU Decomposition: A = LU 

|n 3 + 0(n 2 ) 

Solve: AX = Y for X with A, X , Y € R nxn 

|n 3 + 0(n 2 ) 

Eigenvector Decomposition: A = T AT -1 

8n 3 + 0(n 2 ) 

Schur- vector Decomposition: A = USU T 

8n 3 + 0(n 2 ) 

Solve Lyapunov eqn. PA + A T P = Q A, P,Q (z R nxn 

12n 3 + 0(n 2 ) 


3.5.1 Eigenvector Method 

The procedure we refer to as the eigenvector method is based on the Mcfarlane-Potter- 
Fath algorithm [39]. The procedure for solving (3.209) is. (supressing the dependence 
on x for notational convenience): 
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1. Form the Hamiltonian Matrix: 


H = 


A -BR~ x B t 

-Q -A T 


(3.215) 


The number of operations in this step is essentially the work involved in forming 

BR~ l B T . 


2. Do an eigenvector de-composition of Ti. If {A,B,C} is stabilizable and de- 
tectable and the eigenvalues of H are distinct and we appropriately reorder the 
eigenvectors we get: 


HT = T A 


1 

IN 

i 


1 

1 

> 

o 

i 

i 

to 

t 


0 A 


(3.216) 

(3.217) 


where T is a matrix containing the eigenvectors of H and — A is an n x n diagonal 
matrix containing the stable eigenvalues of "H. 

3. The positive definite solution to the Ricatti equation is then given by: 


-l 


p = t 21 t- 


(3.218) 


4. Finally we compute the feedback gain: 

A'(x) = R~ l B T P (3.219) 

The following table sumarizes the operation counts for this method, where we have 
assumed that B 6 R nxn : 
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Step 

Operations count 

Form H: 


Form: Rr x B T 

i n 2 3 

Multiply: B * (R~ l B T ) 

n 3 

Eigenvector decomposition of H 

8(2n) 3 

Form P = V 2 1 V {[ 1 

±n 3 

3 1 

Compute feedback gain K\ 


Multiply: -(R~ 1 B t )*P 

n 3 

Note: R~ l B T is available from step 1 


TOTAL 

CO 

£ 

cr> 

CO 

?? 


3.5.2 Schur-vector Method 

The procedure we refer to as the Schur-vector method was developed by Laub [38] 
and differs from the Eigenvector method in that it uses the numerically more reli- 
able Schur-vector decomposition [18] in stead of an Eigenvector decompostion. The 
procedure for solving (3.209) is, (supressing the dependence on x for notational con- 
venience): 


1. Form the Hamiltonian Matrix: 


H = 


A -BR~ l B T 
-Q -A T 


(3.220) 


2. Do a Schur-vector de-composition of H. If {A,B,C} is stabilizable and de- 

tectable we get: 


HU = US 


r 

£ 

H 

tS* 

i 


1 

Co 

Co 

to 

i 

£ 

*o 

1 


1 

(N 

(N 

to 

O 

i 


(3.221) 

(3.222) 
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where S is an upper triangular matrix with the eigenvalues of H appearing on 
its diagonal. It is furthermore assumed that the Schur- vectors have been re- 
ordered so that the stable eigenvalues of H appear on the diagonal of 5 X1 while 
the unstable eigenvalues appear on the diagonal of S 2 2 - 

3. The positive definite solution to the Ricatti equation is then given by: 

p = U 2X U; x l (3.223) 

4. Finally we compute the feedback gain: 

A'(x) = —R~ l B T P (3.224) 

The following table sumarizes the operation counts for this method: 


Step 

Operations count 

Form Pt: 


Form: R~ l B T 

3 n 

Multiply: B*(R~ 1 B T ) 

n 3 

Schur-vector decomposition of H 

8(2n) 3 

Form P = S'nSu 1 

4 n 3 

3 U 

Compute feedback gain K : 


Multiply: -(R- l B T ) % P 

R 3 

Note: R~ l B t is available from step 1 


TOTAL 

69n 3 


3.5.3 Kleinman’s Method 

Kleinman’s method [29] is an iterative procudure for solving the Ricatti equation. It 
uses the following algorithm assuming that a starting guess of P, say P 0 is available: 

For i = 0, 1 . . . 


99 




1. Form: K, = —R 1 B T P, 


2. Form: A, = A — BK l 

3. Solve the Matrix Lyapunov equation: 

Aj P l+ 1 + P i+ 1 A { = -{Kj RK i + Q) (3.225) 

for P i+ i 

4. If = P, we have a solution to (3.209). Otherwise iterate this process by 

setting Pi P{+i and starting again at step 1. Note that at the completion we 
already have the feedback gain available from step 1. 

This iterative process should be natural to implement for our problem since we can 
assume that we have a previous value of P(x) computed for a nearby state x, which 
we can use to start the iterations. The following table gives an operation count for 
one iteration: 


Step 

Operations count 

Form 


Multiply: B T * P, 

Ti ^ 
It 

Form: R~ 1 B T P t = K t 

I " 3 

Form Ak\ 


Multiply: B * 

n 3 

Subtract: A — ( BK { ) 

n 2 

Form Kj RK{ + Q 


(Note: KjRKi = {P i B)*{R- 1 B T P i ), so 


use factors from first step) 


Multiply: PiB *(R~ l B T Pi) 

n 3 

Add: KjRKi + Q 

n 2 

Solve the Lyapunov Equation: 

12n 3 

TOTAL for k iterates 

~ k * 16n 3 
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3.5.4 Collar and Jahn’s Method 


The method we refer to as Collar and Jahn's method, is based on the Mcfarlane- 
Potter-Fath algorithm, but uses an iterative method to solve for the eigenstructure 
of the Hamiltonian matrix, in stead of using the more standard Hessenberg-QR- 
method (see e.g. [18]). The steps for Collar and Jahn’s method are the same as for 
the eigenvector method outlined above except that the eigenstructure of H is found 
as follows [11]: 


1. Assume that we have an iterate for the eigenvectors and eigenvalues of H. say 
V k and D k ■ Find A*, and D k+i from: 


Dk+\ 

A k 


diag V k 1 HVk 
V k x HV k - D k+1 


(3.226) 

(3.227) 


i.e. D k+ i consists of the diagonal elements of V k x HV k , and A k consists of the 
off-diagonal elements of the same matrix. 

2. Solve for E k from: 

E k D k+1 - D k+1 E k = A k (3.228) 

Since D k+ j is diagonal, the solution is simply given by: 


[E k ] tJ 


[Afc]»j 

[Dk+i}u + [D k+ i ]jj 


3. Form the updated estimate of V\ 


(3.229) 


V k+l = V k (I + E k ) 


(3.230) 


4. Repeat the process until A*, is small enough. 


The operation count for one iterate is given by: 
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Step 

Operations count 

Form 7i: 


Form: R~ l B T 

i n 3 

3 U 

Multiply B * (R~ x B T ) 

n 3 

Update Eigenvectors of 7i 


Form: V r *T 1 7 i. 

|(2n 3 ) 

Multiply: (V^~ 1 'H) * V* 

(2n) 3 

Solve for A k 

(2n) 2 

Form: 14+1 = Vfc * (I + Ek) 

(2n) 3 

Form P = V^T 1 

- XI 3 
3 

Compute the feedback gain K: 


Multiply — ( R~ l B T ) * P 

n 3 

TOTAL for k iterates 

k * 31n 3 


3.5.5 Recommendation 

Based on the operation counts we have found above we generally prefer the Schur- 
vector method for on-line computation of the feedback gains, since: 

• it is competitive in terms of computational load (assuming the iterative methods 
require 3 to 5 iterations per solution) 

• it is numerically more reliable than the eigenvector method 

• it has a more predictable execution time — it does not have the open ended 
iterative nature of the other two methods. 
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3.6 Simulation Results 


In this section we will apply the continuous Ricatti design method to a nonlinear 
system. We will see that the method generally works better that designs based on 
the linearized dynamics but that at some point it becomes unstable. An approach 
to improve the stability properties will be discussed in Chapter 4 where we will use 
an optimal control approach which utilizes knowledge about the future behavior of 
A(x),tf(x). 

3.6.1 System Dynamics 

The system we use for the simulations will have dynamics related to that of a one-link 
flexible manipulator. We will change the dynamics of a one link arm by simplifying 
the equations of motion on the one hand, but at the same time introduce parameters 
which can be used to accentuate the main nonlinearity in the problem. 

To derive the equations of motion for the one link flexible arm of Figure 3-11 we use 
the assumed modes approach discussed in Appendix A. For this example we assumed 
a single sinusoid mode shape so that the displacement of each “slice” of the arm 
relative to the line joining the root and the tip (see Figure 3-11) is given by: 

«(C>9) = sin (3.231) 

where: 

L = length of the link (3.232) 

£ = “location” of the slice (3.233) 

In what follows we shall refer to 9 in Figure 3-11 as the gross motion angle and to q 
as the flex mode deflection. 
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t 



Figure 3-11: One Link Arm 

The resulting equations of motion can be put in quasilinear form. We get: 

x = A(x)x + B(x)u (3.234) 

with: 

x T = 6 q 6 q 

0 0 10 
0 0 0 1 
0 032 O33 0 

0 042 043 0 



(3.235) 


(3.236) 
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0 


B(x) 


0 


- 8.905 

q q.^|q 4.636x10 3 +'29.61<7 2 
~t(q) 


(3.237) 


where 


032 


033 


042 


043 

7(?) 


(-0.218 x 10 7 + 18. 85(740. 20 2 )) 

7 ( 9 ) 

-18.85<?(3.142g) 

7(7) 

-8.5 x 10 4 (3.6 x 10 2 + 3q 2 ) + 1.733 x 10 2 (1.132fl 2 ) 

7 (<?) 

-1.733 x 10 2 g(1.61 x 10 2 <? - 9.425flg 2 ) 

7(9) 

1.157 x 10 s - 2.961 x 10 2 q 2 


(3.238) 

(3.239) 


(3.240) 


(3.241) 

(3.242) 


and u is the control input. We used parameters similar to the space shuttle R.M.S., 
viz: 


L (length) = 

13.44m 

(3.243) 

pA (density per unit length) = 

55.16 kg/m 

(3.244) 

El (stiffness) = 

1.38 x 10 1 1 *2.08 x 10 ~ 5 Nm 2 

(3.245) 


If we neglect terms like q 2 , qq , q 3 9 and round some of the constants A(x), B(x) simplify 
too: 


A(x) 


0 0 10 

0 0 0 1 

0 -an + a u 0 2 0 0 

0 — a 2 i + ct-nQ 1 0 0 


(3.246) 
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0 


B(x) 


0 

-Pi 

02 


3.24 




where for the nominal case: 


Oil 

= 20 

(3.248) 

Ot 12 

= 0.12 

(3.249) 

Q 2 l 

= -267 

(3.250) 

<222 

= 1.7 

(3.251) 

Pi 

= 8 x 10“ 5 

(3.252) 

02 

= 1.4 x 10 -3 

(3.253) 


To find a feedback law we used the approach discussed in Section 3.2.3. Our goal was 
to cause the gross motion angle 6 and angular rate 6 to remain similar to that of a 
reference model while at the same time keeping the flex mode deflection small. We 
used the cost function and augmented dynamics of equations (3.47), (3.48), (3.53). 

The reference model we used is: 


where: 


= 


A d x d 
0 


— UJl 


1 

-2C^'n 



l 

Qa 

ft* 

i 


. ^ . 


(3.254) 


(3.255) 


io n = 10 


(3.256) 


C = 1 


(3.257) 
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For this example (which is intended to be stressful to the continuous Ricatti design 
method) we used the following parameters in the cost function (see equation 3.53): 



e r = 

(o-e d ) (o-e d ) 

(3.258) 


Z T = 

6 q 6 q 6 d 6 d 

(3.259) 

and: 




Qy 

= diag 

1 x 10 7 lx 10 5 

(3.260) 

Q , 

= diag 

r 

0 1 X 10 5 0 0 0 0 

(3.261) 

p 

= 1 x 10" 4 

(3.262) 

Case 1 




For the first case we used the cost 

function etc. as specified above and 

accentu- 

ated the nonlinearity of the system by using the following values for a tJ 

in equa- 

tions (3.246), (3.247): 


= 20 

(3.263) 


<*12 

= 0.12 

(3.264) 


<*21 

= -26 

(3.265) 


<*22 

= 50 

(3.266) 


,3i 

= 8 x 10" 5 

(3.267) 


$2 

= 1.4 x 10 -3 

(3.268) 

The desired response is a 

90° re-orientation of the “arm” in about 0.5s which is a 

strenuous manever for an 

“arm” of this size. Figure 3-12 shows the stable 

response 


of the closed loop system when using the continuous Ricatti design feedback, while 
Figure 3-13 shows that the closed loop response for the feedback law based on the 
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linearized dynamics is unstable. 


Case 2 

For this case we further accentuated the nonlinearity by using: 


On = 20 

(3.269) 

012 — 0.12 

(3.270) 

a2i — -26 

(3.271) 

022 = 170 

(3.272) 

\r> 

1 

O 

X 

oo 

II 

(3.273) 

02 = 1.4 x 10" 3 

(3.274) 


At this point we see from Figure 3-14 that the continuous Ricatti design feedback 
law results in unstable behavior. We will revisit this case in Chapter 4 where we will 
discuss an approach to improve the stability properties by using an optimal control 
approach which utilizes knowledge about the future behavior of A(x) and B(x). 


3.7 Summary 

In this Chapter we explored a control design method that exploits the quasilinear 
form of dynamics for nonlinear systems. The method can be applied to a large class 
of nonlinear systems and provides the designer with convenient design parameters 
to modify the closed loop system behavior. It has been our experience (see also 
Chapter 5 and [62] for further examples) that the method resulted in well behaved 
closed loop systems if we appropriately adjusted the design parameters. However only 
local stability can be guaranteed. To analyze the “non-local” stability properties we 
showed the connection between the continuous Ricatti design method and receding 
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Figure 3-13: Unstable Response — Linearized Design 
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horizon control. We also showed how Zubov’s method could be used to numerically 
determine the stability domain of the closed loop system. 

In the next chapter we will explore an optimal control methodology that further 
exploits the quasilinear form of the dynamics and can be used to improve closed loop 
stability. 
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Chapter 4 


Controllers with Look-Ahead 


4.1 Introduction 

In this chapter we further explore and exploit the similarity between linear time- 
varying systems and nonlinear systems in quasilinear form. The main difference 
between the methods presented here and the continuous Ricatti design method of 
the previous chapter is that the new control methods take into account “knowledge" 
about the future behavior of A(x) and Z?(x). Such controllers will be refered to as 
“controllers with look-ahead”. 

Sections 4.2, 4.3 and 4.4 examine control strategies for general (not necessarily feed- 
back linearizable) nonlinear systems which will result in stable closed loop systems. 
The methods are all rooted in optimal control theory and require one to solve opti- 
mization problems of the form: 

Given a system: 

X = f(x, u, t), x(t 0 ) = X 0 


y - h(x) 
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(4.1) 

(4.2) 



Find the control input u that will minimize the cost function: 


J = <p(x(t /)) + / £(x,u,r)dr, £(x. u,r)>0 

Jt 0 


(4.3) 


Note that since C is explicitly a function of time, the optimization problem as defined 
above is sufficiently general to include tracking problems. For example assume that 
the plant output y is required to track some reference trajectory, say y d (t), then we 
can use: 

£(x,u,<) = C(x,u,y d (t)) (4.4) 

Once the general theory for the methods is in place, we develop an approach to 
(approximately) solve the required optimal control problems by using the similarity 
between quasilinear systems and linear time varying systems (see Section 4 . 5 ). 

Assumption 4.1.1 We will assume that our systems are sufficiently well behaved so 
that if we apply the optimal control for the problem defined by equations (4-1), (4-2), 
(4-3), then £(x, u,r) will remain bounded over the interval [t 0 , ^ /] and t p(x(tf)) will 
also be bounded. 

Assumption 4.1.1 is needed to ensure that £(x, u, £) does not, for instance, behave 
2 

like ( 7 ^) 3 which becomes unbounded at t = a but remains integrable across the 
singularity. Note that Assumption 4.1.1 does not in itself guarantee that the closed 
loop system will be stable, since some of the system states may well be unstable but 
not affect C. 
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actual trajectory 
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Figure 4-1: Time Intervals for Short Term Tracker 

4.2 Short Term Optimal Tracking 

This section examines an approach to control nonlinear systems by solving a sequence 
of short term optimization problems. The goal of the method is to cause the output 
of the system (4.1), (4. 2) to track some desired trajectory, say y d (t) and at the same 
time ensure that all the state variables remain bounded (see also Section 3.2.3). It is 
furthermore assumed that at any time t the reference trajectory is onlv known over 
a limited interval (horizon), say [t,t + T], 

To achieve the tracking/regulation goal, given the limited knowledge of the desired 
trajectory, the control law is found by solving a sequence of optimization prob- 
lems each defined for a time interval [tiJ, + T], i = 0,1... (see also Figure 4-1). 
Lemma 4.2.1 shows that under certain conditions this approach will lead to stable 
closed loop systems. Before showing the desired stability result we establish the 
following well known fact: 
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Fact 4.2.1 Assume that the nonlinear system: 

x = f(x,u,t) (4.5) 

is completely controllable. Consider the cost function: 

J = px T (tf)Hx(tf) + f £(x,u,r)dr (4.6) 

Jt 0 

where £(x,u,r) > 0 and H is positive semi- definite. Let u M [to,^/] be the control time 
history which minimizes J for a specific value of p, and let x^[fo, tf\ be the associated 

state trajectory. Then if we apply the optimal input u^, to the system, !/) 

will remain finite. Furthermore as p oc we will have x^(t j)Hxfit f) — * 0. 


Proof: 

Since (4.5) is completely controllable there exists some control history, say u[to, tf\, 
which drives the system from its initial condition to a state where x T (t f)Hx(tj) = 0. 
Let x be the associated state trajectory and 


J = 



(4.7) 

(4.8) 


be the cost associated with using u as the control. Since u M is optimal it follows that: 


ft f 

0 <px T Hx^, + / £(x M , u m . r)dr < J 

Jt n 


(4.9) 


Since £(x, u, r) > 0 we then have: 


P*lH x M 

tf 

< 

j 


i e 

< 

j 

i i 


(4.10) 

(4.11) 


Since J is finite x^Hx^ will be finite. Furthermore since J does not depend on p we 
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can make the left hand side of (4.11) arbitrarily small by making y large enough. 


Q.E.D. 


With this fact in hand it can be shown that under certain conditions the strategy 
outlined above will result in stable closed loop systems (in a bounded input-bounded 
output sense). In order to simplify the statement of the proof of stability for this 
method we assume that a state transformation has been applied to the system so 
that the state vector explicitly contains the plant output, i.e: 


x = 


y 


X r 


(4.12) 


Lemma 4.2.1 Consider the nonlinear system (4-5). Assume that the system is 
reachable for any time interval > T r , and that the reference trajectory y d is bounded. 
Divide the time axis into a sequence of time intervals [f,, t z -f T],i = 0,1,... with 
T > T r ( see e.g. Figure 4-1) , and consider a sequence of optimization problems, one 
for each time interval, where the cost function for each interval is: 


J, 


y ( e T H y e + x t t H x x r 


ft f x 

) +/ ‘(e^^e + x r Q x x r + u r f?u)dr 

'/ Jt , 


(4.13) 


where e = y — y d and Q y , Q x , H y , H x are all positive definite. Let u,[ft, t t + T] be the 
control input which is optimal for the i’th problem. If Assumption 4 . 1.1 holds, and 
p is large enough, and we then apply u.ff,, t, + T] on each interval, the closed loop 
system will be stable. 


Proof: 

Consider an arbitrary interval [t,,t, + T}. Assume that the state of the system at 
the beginning of the interval, i.e. x(C), is bounded. Then Fact 4.2.1 implies that, if 
u[C,C + 71 is applied over the interval i /], then e T H y e + xjH x x r will be arbitrarily 
small, provided that y is large enough. Since H y and H x are positive definite, this 
implies that e(tj) and x r (tf) will be finite, and subsequently since is bounded, 
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x(tf) must be bounded. Furthermore since Q y ,Q x are positive definite, Assump- 
tion 4.1.1 ensures that x will remain bounded. At this point we have established 
that: if the system state is bounded at the beginning of an interval, it will remain 
bounded throughout the interval up to the beginning of the next interval. Assuming 
that the state of the system at t 0 is bounded, it follows that the system response is 
bounded for all intervals. 

Q.E.D. 


The result above shows that by solving a sequence of short term optimal control 
problems we can obtain stable closed loop systems. Unfortunately, in practice, this 
approach has limited value since the control history generally is discontinuous at 
the boundary of each interval and the control law is thus likely to excite unmodeled 
dynamics which may lead to instability. (See also Example 4.6.1). 


4.3 Receding Horizon Control 

The previous section examined a stable control law which required that a sequence of 
short term optimization problems had to be solved. The method had the drawback 
that it would lead to discontinuous control histories. To remedy the problem of the 
discontinuous control we examine another control methodology, i.e. receding horizon 
control (see e.g. [36], [41]). 

The receding horizon control strategy is: 

1. At the current state, say x(<) solve an optimization problem of the type defined 
by equations (4.1), (4.2), (4.3), over a finite horizon [ t,t + T]. 

2. If u[f , t + T] is the solution to the optimization problem over the full interval, 
apply only the initial control, i.e u(f), at the initial state. 
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3. Repeat the process starting at step 1. 


Fundamental results regarding the stability of this control law are available. Kwon 
and Pearson [36] show that for linear time varying systems: 

x = A{t)x + B(t)u (4.14) 


with a cost function of the form: 

ft+T 

J = Jt (x T Qx + u T /?u)dr (4.15) 

and subject to the terminal constraint that x(t + T) = 0 the closed loop system 
will remain stable. Mayne and Michalska [41] expand this result to nonlinear sys- 
tems. They show that if the receding horizon control law is applied to systems of the 
form (4.1) with cost function: 

J = ^ (x t Qx + u r i?u)dr Q,R> 0 (416) 

and subject to the terminal constraint that x(t + T) = 0 the closed loop system 
will be stable. Note that in both the linear and nonlinear system cases the above- 
mentioned stability results cannot be applied as is to tracking problems, since for 
tracking problems y^t + T) / 0 in general. Stability results for tracking problems 
appear in [34] for linear time invariant discrete time systems: 

x(k + 1) = Ax(k) + Bu(k) (4.17) 

y = Cx(k) (4.18) 

and cost functions of the form: 

i i j=fc+iv 

J=-(e T (k + N)He(k + N)) + - £ (e T U)Qe(j) + uU) T Ru(j)) (4.19) 

“ j = k 
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where: 


e (j) = y(i) — Yd(j) (4.20) 

and it is assumed that the reference trajectory yd(j) is available over the horizon 

[k,k + N\. 

In what follows we extend the stability results mentioned above to the nonlinear 
receding horizon tracking problem (he. where x(f + T) ^ 0). In the process we also 
develop an alternate method of proof for the nonlinear receding horizon "regulator” 
stability result by Mayne and Michalska (i.e. where x(t + T) = 0). Our results depend 
on the following Lemma: 


Lemma 4.3.1 Consider the nonlinear system (4-1) and cost function (4-3). Assume 
that the system is reachable and allow for the possibility of a terminal constraint of 
the form xj>(x(tf)) = 0. Let u[to,t j] } and x[t 0 ,tf] denote the optimal control and state 
trajectories respectively and let: 

V{x,t 0 ,tf)= [ / (£(x(r),u(r),r)d< (4.21) 

Jto 


denote the actual cost incurred if we apply the optimal control u over the interval 
[to,tf] starting from state x at time i 0 * Then: 


-V(x,t,t + T) 


-£(x, u, () + 


x=f(x,u,t) 


dVjx, tp ,tf) 

dtf 


t f =t+T 


(4.22) 


Proof: 

We have: 


dK(x, t, t + T) 

At 


dV{x,t 0 ,t f ) 

dx 


dx dV(x,t 0 ,tf) 
t 0 =t d t dto 

t f = t + T 


dv(x,t 0 ,tf) 

dt } 


<o = t 
t/ = f+r 


t 0 = t 
t/=t+r 


(4.23) 
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dV(x,t Q ,t f ) 


+ 


dx 

dV(x,t Q , tf) 


t 0 =t 
t/ = t+r 


f(x(o,u(t),o+ — 

Ot Q 


<o = ' 
i/ = l + T 


dt 


f 


t 0 =t 

t f =t+T 


(4.24) 


From the HJB-equation we know that: 


dV{x,t 0 ,tf) „ dV(xJ 0 Jf)„, 

= C(x, u(< 0 ), t 0 ) + c)x ~ f (x^ u(/ 0 ), io) 


(4.25) 


Substituting (4. 25), into (4.24) then gives us equation (4.22) as desired. 


j Q.E.D . 


We can find an expression for 


as follows: 


dV(x,t 0 ,tf) 

dtj 


(4.26) 


Fact 4.3.1 Consider an optimization problem where we have a system as in (4.1). a 
cost function: 

J = ¥>(x( */),</) + / 1 £(x,u,r)dr (4.27) 

Jto 

and a terminal constraint ^(x(t f ),t f ) = 0. Let V{x,t 0 ,t f ) again be the actual cost 
incurred if we apply the optimal control. Then: 


dV(x,t 0 ,tf) ( d$(x,t f ) d${x,t f ) N 

W, + ~mr~ + £(x - t '\ 


I '/ 
*(*/) 


where: 


*(x,tf) = ^(x(</),i/) + <p(x(tf),tj) 
and v is a Lagrange multiplier. 


(4.28) 


(4.29) 


Proof: 

Let u[fo,ty] and x[£o*^/] be the optimal control and state trajectory associated with 
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the optimization problem defined above, i.e. over [to,tj]. Define the augmented cost 
function: 

J a = f v T jj,(x(t f ),t f ) + v(x.(tf))+ f’ (£(x,u,r) + A T (-x + f(x,u,r))}#.30) 

J t o 

= f $(x(f/),f/) + [ S (?f(x,u,r) - A T x) dr (4.31) 

Jto 

where u, A are Lagrange multipliers and we have defined: 

7-{(x,u,r) = f £(x, u, r) + A r (f (x, u, r )) (4.32) 

d = { u T i>(x(t } )) + 0(x(tj)) (4.33) 

Now consider the perturbation in J a if we perturb: 


• the terminal time to: tj + At 

• the control history to: u(f) + Su(t). Assume that this perturbation in the 

control history results in a perturbation of the state trajectory to: x(f) + 6x(t). 
Furthermore assume that the perturbed state trajectory satisfies the terminal 
constraint: i/>(x T | ^ ^ = ® and that the initial condition x(to) remains 

unchanged. 


The perturbation in the augmented cost, after integration by parts of the A term, is 
then given by (see e.g. [8]): 



(4.34) 


By appropriate choice of the Lagrange multiplier A, and the necessary conditions for 
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u to be optimal, the only remaining terms are: 


' M(x ’ < h i + ^h + £(x , u , /j) ' 


dx 


dt< 


At 


*u/> 


(4.35) 


So that: 


dV{x, to, if) 
dtf 


= lim 


A J a 


&t-~o At 

'd$(x,tj) . d$(x,t f ) 


dx 


-X + 


dt i 


+ £(x, u, t f ) 


(4-36) 

(4.37) 


z t 

*(t f ) 


Remarks: 


Q.E.D. 


• Given assumption 4.1.1, will be bounded. 

• We expect that when tf — to becomes large, the change in cost A J a and thus 
|^- will become relatively small. In equation (4.28) this should occur through 
u changing its value. 

• We can recover Michalska and Mayne’s result by noting that their result is 
derived for the case where: 

- The system is time invariant with a single equilibrium point at the origin, 
viz: 


x = f(x,u) (4.38) 

0 = f(0,0) (4.39) 

- There is no explicit terminal penalty, i.e.: ^(x(tf) = 0 

— The terminal constraint is: 


y>(x, tf) = x(t f ) --- o 


(4.40) 
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- The integrand of the cost function is: 


C(x , u) = x T Qx + u T Ru 


(4.41) 


where Q, R are both positive definite. 
Under these conditions we see that: 


- V (x, t, t + T) = V (x) i.e. the actual cost to drive the system to the origin 
starting at state x at time t , depends only on x. 

- V(x) is a radially unbounded positive definite function and can serve as a 
candidate Lyapunov function. 

- Because of the terminal constraint x(tj) = 0, and form of the cost function 
we have for case where x = f(x) + G(x)u that u (tj) = 0 and x(f/) = 0 
and thus: 


dVjxJo ,_tf) 

dtf 


d$(x,t f ) . d$(x,t f ) 


dx 


-x + 


dt, 


+ C{x,u,t f ) 


Vf 

M*f) 


(4.42) 


= 0 


(4.43) 


— Stability then follows by noting that under these conditions: 

—V(x,t,t + T) = ~{x t Qx + u t Ru) 

< 0 Vx # 0 


(4.44) 

(4.45) 


We can now proceed with the stability result for the nonlinear receding horizon track- 
ing problem. Our goal will be to show that if we apply the receding horizon control 
law, then “on average” the full state will remain bounded. We take this approach 
because in general it is not possible to cause the system (4.1) to perfectly track the 
reference trajectory and at the same time ensure that all the system state variables 
remain bounded (see e.g. [59]). This means that it generally will not be possible to 
show that the error variables, say e = y-y d , and the “rest of the state”, say x r , tend 


124 



to 0 as t 


oo. 


Lemma 4.3.2 Assume that we have a system of the form (f.l), and that we have 
chosen a cost function appropriate for a tracking problem, say: 

rt+T 

J = J £(x,u,<)dr (4.46) 

rt-\-T 

= J (e T Q y e + x r Q x x r + u T Ru)dr (4.47) 

where e(r) = y(r) — y d(r), and Q y ,Q x ,R are all positive definite and we have a 
terminal constraint such that (e T Q y e + x r Q x x T ) tf = 0. Then if Assumption 4 . 1.1 
holds we will have: 


1 c v 

lim - / ( e T Q y e + x r Q x x r + u r Ru)dr < k 
'J— 00 r > Jo ’ ~ 


(4.48) 


for some k 


Proof: 

Let V (x, t,t + T) be the cost incurred if we apply the optimal control over the interval 
[M+T] starting at the state x. Then if we apply the receding horizon control strategy, 
we know from Lemma 4.3.1 that: 


d_ 

d t 


V(x,t,t + T) = — £(x, u, t) + 


dV(x,t 0 ,tf) 

dtf 


i 0 =t 

t/=t+r 


(4.49) 


Since assumption 4.1.1 holds we know that < k for some n > 0, so that 

equation (4.49) becomes: 


—V(x,t,t + T ) < — £(x, u, t) + k 


(4.50) 


Integrating both sides of (4.50) gives: 


[V dV p p 

Jo d? dr - ~ Jo £ (x,u,T)dr+y o xdr (4.51) 
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V'(x(r?), 7/, ?7 + T) - F(x(0),0,T) < - / (e r Q„e + xjQ x x r + u r #u)d7 

V 0 


+ / /tdr 

Jo 


(4.52) 


Now by our choice of £ we know that V(x,t,t + T) > 0 Vx, t. Thus: 

— V(x(0), 0, T) < - / (e T Q y e + x*Q x x r + u T Ru)dr 

Vo 

+ f Acdr (4.53) 

=>• r (e T Q y e + xjQ x x r + u T i?u)dr < V(x(0),0,T) + / fcdr (4.54) 

Vo 

1 rn T r . .■ /V(x(0),0.T) 

=>\im- (e T Q y e + x J r Q x x T + u r i?u)dr < lim 

r/ Jo \ V 

+ — f Kdr J (4.55) 

V Jo ) 

=*> lim - / (e r Q v e + x^Q x x r + u T /?u)dr < k 

v^oo n j o 


(4.56) 


j Q.E.D . 


Lemma 4.3.2 provides us with an “integral boundedness” result. We see that “on 
average” e T Q y e + x*Q x x r has to be smaller than «. Since y d is assumed to be 
bounded, this then implies that x will “on average” have to be bounded. 


4.4 Long Term Optimal Control 

The receding horizon controllers discussed in the previous section provide continuous 
control laws but require much computation: they require that one repeatedly solve 
an optimal control problem. We have found that the more traditional long term 
optimal controllers, as defined by equations (4.1), (4.2), (4.3) generally are more 
useful. These long term controllers are typically used to transfer a system from some 
initial state to a new equilibrium point. Such will be the case if we want a flexible 
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robotic manipulator to perform a large scale maneuver such as moving a payload 
from one location to another. Note that the long-term optimal controller requires 
one to solve an optimal control problem only once for each maneuver. The receding- 
horizon control law on the other hand, would in principle require an infinite number 
of solutions to the optimal control problem for the single maneuver. 

The stability of long-term optimal control laws are guaranteed bv the following re- 
sult [2], [45] which also shows that the optimal control laws are robust, even for the 
case where we are dealing with nonlinear systems. 

Theorem 4.4.1 Consider the system: 

x = f(x) + G(x)u (4.57) 


and take as a cost function: 


J = / ( m (x) + u r Ru)dr (4.58) 

J to 

where it is assumed that R > 0 and m(x) is a radially unbounded globally positive 
definite function . Assume that: 

• The system (4.5 7) is completely controllable. 

• The functions f (x), G(x), m(x) are sufficiently well behaved such that an optimal 
control exists and that the optimal control satisfies the HJB-eyuation. 

• The free system: 


f(x) 

(4.59) 

m(x) 

(4.60) 


is completely observable . 
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Then if we apply the optimal control input, say u, the closed loop system will be 
globally asymptotically stable. Furthermore if we apply the (incorrect) control input: 

u = u + k(t) u (4-61) 

the closed loop system will still be asymptotically stable provided that: 

^ < k(t) < oc (4.62) 

4.5 Look-Ahead Control Using the Quasi-Linear 
Form 

4.5.1 Introduction 

The control methods discussed in Sections 4.2, 4.3 and 4.4 all require us to solve 
optimal control problems for a nonlinear system. In this section we construct (and 
use) approximate solutions to the nonlinear optimization problems by exploiting the 
quasilinear form for the dynamics. We do this by performing the following steps: 

1. Convert our system dynamics to quasilinear form, say: 

x = f(x) + G(x)u (4.63) 

= A(x)x + i?(x)u (4.64) 

2. Estimate the trajectory that the system will follow, for example x(f) « x(t). 

3. Treat the quasilinear system as if it were a linear time varying system along the 
estimated trajectory, viz: 


x = A(x)x + B(x) u 


(4.65) 
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^ A(x(f))x + B(x(t))u 
= A(t)x + B(t)u 


(4.66) 

(4.67) 


4. Solve the (much easier) optimal control problem for the linear time varying 
system. 

5. Apply the feedback law found from the linear time varying problem to the 
nonlinear system. 


Comparison with standard linearization 

Note that the approach outlined above is different from the usual method of linearizing 
a system along a reference trajectory, and then finding a controller for the linearized 
system. To see this recall that the process of linearizing a system around a reference 
trajectory involves the following steps: 

• Assume that we have a nonlinear system: 

x = f(x) + G(x)u (4.68) 

• Let u n be a nominal input which causes a nominal trajectory x n , i.e.: 

x n = f(x n ) + G(x„)u n (4.69) 

• Now consider applying an input u n + <5u to (4.68) and write the resulting tra- 
jectory as x n + fo. Equation (4.68) becomes: 

x n + Sx = f(x n + 6x) + G(x„ + <5x)(u + 6u) (4.70) 
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• If Sx remains small, we can use a Taylor series expansion for the nonlinear 
functions. If we retain only first order terms we get: 


x n + Sx 


=> Sx 


f (x„) + 


df' 


j)xj 
+G(x n )u n + 


Sx 


Xn 


d_ 

dx 


(G(x)u„) 


<^x + G(x n )Su 


df 9 

-+^(G(x)u„) 


Sx + G(x n )<5u 


(4.71) 

(4.72) 


• The linear time varying system that we obtain through linearization is given 


by: 


Sx = A(t)<5x+ B(t)Su 


(4.73) 


where: 


4(1) 

B(t) 


def 


M 9 <ru x x 

dk + *; (G(X)U " ) 


= G(x n ) 


(4.74) 

(4.75) 


• We now find a control law for the linearized system, say: 


S u = K(t)Sx 


(4.76) 


• To control the full nonlinear system we then add the perturbation control input 
(Su, to the nominal control, viz: 

u = u n + K(t)Sx (4.77) 

The essential difference between our approach and the standard linearization approach 
is that the latter requires u n to be known a priori, while the approach based on the 
quasilinear form will generate the full u in one step. 
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4.5.2 The Look-ahead Issue 


As indicated earlier the main difference between the methods of this chapter and the 
previous chapter is that we now use knowledge about the future behavior of ,4(x) 
and B(x) in the control laws. To clarify the need for look-ahead let us examine the 
linear time varying quadratic regulator problem where the system dynamics are: 

x = A(t)x + B(t) u (4.78) 

and the cost function is: 

J = x{t f ) T Hx(tj) + f \x t Qx + u T Ru)dr (4.79) 

Jto 

It can be shown [28] that the optimal control input is given by: 

u(t) = -R~ l B(t) T P(t,tf)x (4.80) 

where: 

(4.81) 

(4.82) 

he. to obtain the optimal control input, we have to calculate P(t,tf) by integrating 
the differential equation (4.81) backwards in time from the final time tf, to the current 
time t. This clearly shows that the control at the current time t depends on the future 
behavior of the matrices A{t), B(t). (Note that this dependence of the optimal control 
input on the future behavior of the system is also reflected for the more general case 
by the fact that the HJB-equation has a boundary condition imposed at the terminal 
time, see e.g. Section (3.2.1)). 

Having established the need for look-ahead we have to address a further issue: 


dP(Tjf) 

foT 11 = Q + P(T,tj)A{T) + A T (T)P{T,t } ) 

t }) B( t) R~ l B t (t )P(r, tj) 

= H 
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The fundamental difference between linear time varying systems: 


x = A(t)x + B(t) u (4.83) 

and quasilinear systems 

x = A(x)x + £?(x)u (4.84) 

is that for linear systems the time histories A(t),B(t) are always the same no matter 
what initial condition we use. On the other hand for the quasilinear systems the time 
histories for A(x(f)), B(x(t)) will generally be different for each initial condition — 
we can think of the quasilinear system as an ensemble of linear time varying systems. 

It is this “unpredictability” of A(x),B(x) that in general makes it more difficult to 
obtain controllers for quasilinear systems. To overcome this problem of not knowing 
the future behavior of A(x), B(x) we will: 

• Define an optimal control problem which will cause the system state to follow 
as closely as possible a specified trajectory, say x<*. 

• Assume that this desired trajectory is achievable and use it to estimate the 
behavior of A(x), Z?(x),viz: 

A(x) « A(x<i) (4.85) 

B(x) « B(xd) (4.86) 

In essence it is by defining/solving an appropriate optimization problem that we are 
able to estimate the future behavior of the system and thus use look-ahead. 


4.5.3 Implementing the Look-ahead Controllers 

In this section we further detail the method used to implement the look-ahead control 
laws. We will focus the method we have found most useful, i.e. to (approximately) 
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implement a long term optimal controller (see Section 4.4) which transfers the system 
from an initial state to an equilibrium point. A brief discussion on how we can use 
a similar approach for the short-term tracker and receding horizon tracker problems 
will conclude the section. In order to more easily discuss the methodology we will 
again assume (although it is not required in general) that the state vector explicitly 
contains the output variables, viz: 


x = 


y 


x r 


(4.87) 


The look-ahead control laws we use, have a similar objective to the controllers dis- 
cussed in Section 3.2.3, i.e. we want the plant output variables, say y, to track desired 
values, say y d , and at the same time ensure that the rest of the state variables x r 
remain small. With this goal in mind we construct a quadratic optimization problem 
similar to that discussed in Section 3.2.3, viz: 

Assume that the desired trajectory is generated by an autonomous system (typically 
linear time invariant) : 


Xrf = A d x d 

(4.88) 

y d = C d x d 

(4.89) 


Form the augmented system dynamics: 

z = F(z)z + G(z)u (4.90) 

where: 


Z — 

X 

F = 

O 

G = 

' B(x) ' 




1 

o 

i 


0 
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Find the control input that minimizes the cost function: 


J 


= z T (tf)Hz(tf) + f } ((y - yd) T Qy(y - yd) + xjQxX-r + pu T Ru)dr(4.92) 

JtQ 

= f z T (tf)Hz(tj) + f \z T Qz + pu T Ru)dr (4.93) 


The look-ahead control law is then implemented by finding an approximate solution 
to this problem. We first estimate z (t) by say z (<), for the interval [fo, </] , and then 
integrate the Ricatti ODE 

^ t,)F(i{r))+F r (z(T))P(r, t,)--P(r ,t,)G(i(r))R-'G T (Hr ))P(T,t,) 

OT p 

(4.94) 

backwards in time starting from the terminal condition: 

P(t f ,t f ) = H (4.95) 

and then apply the (approximately optimal) control input: 

u = --R~ 1 G T (z)P(t. t f ) z (4.96) 

P 

to the full nonlinear system (4.90). 

Remarks: 

• By using a cost function of the form (4.93) and letting p —* 0 the system will 
follow the desired trajectory and at the same time keep x r small, if this is at all 
possible. Thus when p is small we will generally be able to better predict the 
actual trajectory. 

• Since yd in equation (4.89) will be a sum of complex exponentials we can gen- 
erate a very wide variety of desired trajectories by appropriately changing the 
dynamics in equation (4.88). 
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• If our goal is to transfer the system from an initial condition to an equilibrium 
state, we have found it expedient to do the following. Use tj^oT where T is the 
largest time constant of the trajectory generator (4.88), and use H = P where 
P is the positive semi-definite solution to the steady state Ricatti equation: 

0 = Q + PF{ 0) + F t {0)P - -PG(0)R- 1 G t (0)P (4.97) 

P 

This can be explained as follows. Transferring the system to a new equilibrium 
condition requires that — ► 0 as t — ► -do. If the controller is capable of 

keeping x r relatively small, we expect that at t & 5 T we will have the full state 
z ~ 0. Once we are at z % 0 we can use a stable controller based on the 
linear approximation of the plant dynamics. As discussed in Chapter 3, this 
approximation near the origin is given by: 

z = F( 0)z + G( 0)u (4.98) 


By using: 

u = —~R~ l G(0) T Pt. (4.99) 

P 

for all t > 5 T we are then simply using the Linear Quadratic Regulator feedback 
law based on the linearized dynamics. Note that the choice of H — P ensures 
that the nonlinear control that was used for t < 57" tends smoothly toward the 
linear control law used for t > 5 T. 

• One method to estimate the behavior of z (t) is to set: 


z — 


Yd 

0 


(4.100) 


since we expect x r « 0 and y % y^. However we have found that for the long 
term optimal controllers the following method provided better results: 
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- At the initial condition, say z 0 , solve the steady state Ricatti equation: 


0 = Q + PqF(zq) + F(z 0 ) t Po PoG(zq)R 1 G t (zq)Pq (4.101) 

P 

— Estimate the future trajectory of z by calculating the response of the linear 
time invariant system: 

4 = (^F(zo) - -G{z 0 )R- l G T (z 0 )p)j z (4.102) 

(Note that this method of estimating z is consistent with the theorem by Koko- 
tovic et al discussed in Section 3.3.3). 

Look-Ahead Controllers for the Short-Term Tracking and Receding Hori- 
zon Problems 

For both the short-term tracking problem (Section 4.2) and the receding horizon 
control laws (Section 4.3) we have to (repeatedly) solve an optimal control problem 
subject to a terminal constraint. For these cases we have found a penalty function 
approach [15] to be useful. The penalty function method typically uses a cost function 
of the form: 

J = nj V>(x) 7 V(x)| + [ * C(x.,u)dr (4.103) 

H f Jto 

and ensures that the terminal constraint V>( x ) = 0 is (approximately) met by making 
fif large. Given this form of the cost function we can then compute the approximate 
optimal control law in the same manner that we used for the long-term look-ahead 
controller as discussed above. 
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4.6 Simulation Results 


In this section we will apply the look-ahead controllers to the one-link arm example 
that was discussed in Chapter 3. We see that by using approximate look-ahead, all 
three methods, i.e.: 

• short-term optimal control (Section 4.2) 

• receding horizon control (Section 4.3) 

• and the long-term look-ahead controller (Section 4.4) 

result in stable responses for the case where the continuous Ricatti Design was un- 
stable. 

In the following examples we will use the system dynamics as in Section 3.6, where 
we had: 


z = F( z)z + G(z)u 


(4.104) 


with state vector: 

z = 6 q 6 q xj 

and dynamics matrices: 


F{ z) 

where: 


0 1 0 
0 0 1 
Oil + Oi2^ 2 0 0 

021 + CH 22 O 2 0 0 



A(x) 0 1 B 

G( z) = 

0 A d J [0 


(4.105) 


(4.106) 


(4.107) 
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0 


£(x) 


Ad 


0 

-d\ 

—02 

0 1 

-u; n ~ 2 


(4.108) 


(4.109) 


We used the values of a tJ , /3,-, u> n and £ which resulted in an unstable continuous 
Ricatti design controller, viz: 


<*11 

= 20 

(4.110) 

<*12 

= 0.12 

(4.111) 

<*21 

= -26 

(4.112) 

<*22 

= 170 

(4.113) 

01 

= 8 x 10~ 5 

(4.114) 

02 

= 1.4 x 10 -3 

(4.115) 


= 10 

(4.116) 

C 

= 1 

(4.117) 


Example 4.6.1 [Short Term Optimal Tracker] 

Figure 4-2 shows the simulation results when we used the look-ahead method to find 
an approximately optimal short term tracking control input. We used a similar cost 
function to Section 3.6, Case 3, except that now we added a terminal penalty so that 
the terminal constraint would (approximately) be met. The cost function for this 
example w'as: 

J = Hf(e T Q y fe + z t Q z jz) + f (e T Q y e + z T Q z z + p 2 u 2 )dr (4.118) 
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where: 


T 


P 

Qyf 


Qzf 


Qy 

Qz 

p 


0.2 

100 


di a g 

diag 

diag 

diag 


1 xlO 7 lx 10 5 j 

0 1 x 10 5 0 1 x 10 3 0 0 

1 x 10 7 lx 10 5 

0 1 x 10 5 0 0 0 0 1 


1 x 10" 4 


(4.119) 

(4.120) 

(4.121) 

(4.122) 

(4.123) 

(4.124) 

(4.125) 


We see from Figure 4-2 that the system’s response reflects the discontinuity in the 
short term tracking control law at the boundaries of the time intervals, i.e. t = 
0.2s, t = 0.4s, t = 0.6s. 


Example 4.6.2 (Receding Horizon Tracker) 

Figure 4-3 shows the simulation results when we used the look-ahead method to find 
an approximation to the receding horizon control law discussed in Section 4.3. We 
used the same cost function etc. as in the previous example except that the terminal 
penalty was changed to: 


II 

o 

o 



(4.126) 

Qyf = diag 

1 x 10 7 0 


(4.127) 

Qzf = diag 

0 1 x 10 5 

0 0 0 0 

(4.128) 




(4.129) 


For this example we used a look-ahead interval of T = 0.2 s which is two times the 
time constant of the desired dynamics. Note that we obtained a smooth closed loop 
response from the system. However the computational cost was large. In principle 
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Figure 4-2: Responses for Short Term Tracker 
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the receding horizon method requires us to repeatedly solve the optimal control prob- 
lem for each state as we move along the trajectory. In practice we re-compute the 
optimal control at short intervals along the trajectory as in this example, where we 
re-computed the optimal control each 0.01.S. 


Example 4.6.3 (Long-Term Look-ahead) 

Figure 4-4 shows the results when we used the long-term look-ahead controller dis- 
cussed in Section 4.5.3. We used the same cost function as in Section 3.6, Case 3. 
To estimate the trajectory we used the method discussed in Section 4.5.3. i.e. we 
solved the steady state Ricatti equation at the initial state, computed the response 
for the linear time invariant system of equation 4.102 and used z (t) as our estimate 
of the system trajectory. We then integrated the Ricatti differential equation (4.94) 
backwards in time. To start the integration process we used P(tf,tj) — P where P 
was the solution the steady state Ricatti equation at the terminal state, i.e. z = 0. 


4.7 Summary 


In this Chapter we examined methods for controlling nonlinear systems that utilized 
"knowledge ,, about the future behavior of the system, viz: 

• Short term optimal control. 

• Receding horizon control. 

• Long-term optimal control. 

All three controllers gave stable system responses for a case where the zero look-ahead 
control law discussed in the previous Chapter was unstable. Of the three methods 
we preferred the long-term optimal control approach, because the short term optimal 
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Responses for Receding Horizon Tracker (Accentuated Nonlinearity) 
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controllers resulted in discontinuous control laws, while the receding horizon control 
laws were undesirable because of the amount of computation required. 
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Chapter 5 


Control of Flexible Manipulators 


5.1 Introduction 

In this chapter we use the methods discussed in Chapters 3 and 4 to develop controllers 
for a complex nonlinear dynamical system: a flexible manipulator similar to the space 
shuttle Remote Manipulator System[20]. Since flexible manipulators are in general 
not feedback linearizable the standard methods for controlling rigid manipulators 
cannot be applied (we assume that the system will have more degrees of freedom 
than actuators). Several researchers (see e.g. [48], [6], [10], [49]) have investigated the 
problem of controlling flexible manipulators but at present no mature methodology 
is available that will guarantee global closed loop asymptotic stability. The methods 
we use to control the system also do not guarantee global closed loop stability, but 
do provide us with a generic approach to deal with the full nonlinear nature of the 
system and have resulted in well behaved closed loop systems. 
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5.2 System Dynamics 


To derive the equations of motion for the flexible manipulator we used the method de- 
scribed in Appendix A, i.e. we used Lagrange’s equations and modeled the flexibility 
using an assumed modes method. We made the following assumptions: 


• The manipulator has two flexible links. 

• The base of the manipulator is stationary in inertial space. 

• The gross motion angles, 0 \, 0 2 measure the relative angle between the links (i.e. 
we use the “relative-co-ordinates” defined in Appendix A). This means that 0 ^ 
measures the angle of the inboard link relative to an inertial reference frame, 
and 0 2 measures the angle of the outboard link relative to a line tangent to the 
tip of the inboard link. 

• The structural flexibility is modeled using an assumed modes approach. Each 
link has one flex-mode with a sinusoidal mode shape, i.e. the deflection of a 
“slice” on a link (see figure 5-1) is given by: 


v{(,q) = q<p( C) (5.1) 

(5.2) 

A consequence of using these mode shapes is that the gross- motion angles 0 \, 0 2 
represent angles measured relative to imaginary lines joining the root and tip 
of each link. Thus by controlling 9 \ and 0 2 we will be directly controlling the 
tip motion. 

• The arm has a payload of 3000/6 attached to the tip of the outboard link. The 
payload is considered to be a point mass, with no rotational inertia. 

• The effects of gravity are not included in the model since it is assumed that the 
manipulator will be operating in space. 
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• Other physical parameters are [27]: 

Stiffness (El) = 

Density/Unit Length (pa) = 
Length of Link \ (L\) = 
Length of Link 2 ( L 2 ) = 


1.378944 x 

■V 

10 11 * 2.081157 X 10 -5 — (5.3) 

m 2 

55.16— 

m 

(5.4) 

6.37m 

(5.5) 

7.1m 

(5.6) 


• The control inputs to the system are torque motors located at the root of each 
link. 


Using the assumptions above, the equations of motion for the system were automat- 
ically generated using the routines described in Appendix A. The resulting system 
dynamics were obtained in quasilinear form, viz: 


x = A(x)x+ B(x) u 


(5.7) 


where the state vector is: 


Ox 


gross motion angle - link 1 

<7i 


flex mode deflection - link 1 

0 2 


gross motion angle - link 2 

92 


flex mode deflection - link 2 

Ox 


gross motion angular rate - link 1 

4\ 


flex mode deflection rate - link 1 

62 


gross motion angular rate - link 2 

4 2 


flex mode deflection rate - link 2 


and the system matrices are given by: 


(5.8) 


A(x) 


0 I 

H' l (q)K H~'(q)C(q,q) 


B(x) 


0 

H~ l ( ;q)T(q) 


(5.9) 
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with: 


#(q) = 

Inertia matrix 

(5.10) 

K = 

Stiffness matrix 

(5.11) 

C(q,q) = 

“Coriolis 1 ' Matrix (see Appendix A) 

(5.12) 


and T(q) is determined by examining the generalized forces [17] acting on the system. 
For the parameters given above, the linearized system dynamics has four poles at 
the origin, and two sets of undamped poles with frequencies % 8Hz and as 12Hz 
respectively. 

The equations of motion above appear to be quite simple, but in fact represent com- 
plex dynamics as can be seen by examining equations (5. 13), and (5.14) which give 
the numerator and denominator of the second-row, third column element of the in- 
verse of the inertia matrix. It is this complexity that makes the routines described 
in Appendix A useful, since they can produce output in symbolic form suitable for 
documentation purposes, as well as directly generate “C” or Fortran code that can 
be, for example, used in simulation models. 


num[i/ 1 ] 23 = — (— 72 7r 3 + 12 7r 5 ) p 2 A Ll + 36 7r°m 2 Z> 2 } q\ x 

+ {(36tt 3 + 12 7r 5 ) p 2 A L\ + ((.36tt 5 - 288 tt 3 ) p \L 2 

+36 pa 7r 5 m 2 ) Zu} <? 21 - j-144p^L 3 7r 2 sin(0 2 - ~j^-) 

t L 1 

+ (^36 7r 4 sin(2 0 2 + 288 it 2 sin(- — 1311. _ 2 0 2 ) 

V V L i L i 

-36 7r 4 sin(^ ii - 20 2 ) - 288 tt 2 sin(2 0 2 - 
-144 L 2 TT 4 p A m 2 sin( ^ ?11 - 2 0 2 )j L\ j q 2 x 
— | ^ — 288 7T cos(0 2 ^r— ) + 71-3 cos($ 2 P 2 ^ 2 

+72 p A m 2 7t 3 cos(-0 2 + -^-)} L\ 
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— {(— 144 7r -h 8 7r 5 — 24 x 3 ) p 2 A L 2 2 + (24x 5 m 2 + 72x 3 m 2 ) p A £ 2 } 

+ ( (^48 x 3 cos(2 0 2 — - — p— ) — 9 x 5 cos(2 0 2 — 11 ) 

l\ -Li L\ 

+ 15 7r 5 — 96 x 3 ) p^L 3 + ^60 x 5 m 2 — 36 x 5 m 2 cos( ‘ — 2 # 2 ) 

+288 m 2 x 3 cos( — p— — 2 # 2 ) — 144 x 3 m 2 
-01 

—288 m 2 x 3 cos(2 d 2 — " ^ n - )^ /3^ Z -2 


— ( —36 mjTr 5 cos( 


2 7T q n 


— 262 ) + 36m 2 x 5 


den[i/ 1 ] 23 = {9/?^I^7 t 4 9 21 + ((6 x 4 - 36 x 2 ) p\L\ + 18^x 4 m 2 1 2 ) ^ 1 } 9n 

+ {(6x 4 — 36 x 2 ) p\L\ + ( (l8 7r 4 — 144 7 r 2 ) p A L 2 + 18 /? 4 X 4 m 2 ) L\ } qh 

— j^— 18x 3 sin( 2 —2 0 2 ) + 18 7r 3 sin(20 2 — - ^ 1 - ) 

— 144 7r sin(2 0 2 — - ^ X1 ) + 144 x sin( 2 — 2 0 2 )^ p^Z, 2 

-72 p 2 A L 2 x 3 m 2 sin( “ 71 ~ ?n -20 2 )}l 3 ? 2 i 

L\ J 

+ {(4 7t 4 + 144 — 48 x 2 ) /)^I 2 — (-72x 2 m 2 + 12x 4 m 2 ) p A L 2 ]• 1 4 
+ | 48 x 2 + 24 x 2 cos(2 0 2 — ) 

—9 x 4 cos(2 0 2 - 2 ^ 9ll )l/2^) p^L 3 

— f— 144 m 2 x 2 cos( 20 2 — ) — 18x 4 m 2 cos( — 20 2 ) 

V L\ L\ 

+30 x 4 m 2 — 72x 2 m 2 + 144 m 2 x 2 cos( 2 — 20 2 )') / 94 Z 2 

Lj / 

+ ^lSm^x 4 — lSm^x 4 cos ~ 20 2 )^ ^ L 2 1 I 3 (6.14) 
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5.3 Simulation Results 


To control the flexible manipulator described above, we used the two methods we 
found most useful, i.e. the continuous Ricatti design method described in Chapter 3 
and the long-term look-ahead controller described in Sections 4.4 and 4.5. In both 
cases we included terms to obtain integral control. 


Basic Optimization Problem 

The basic structure of the controller is derived from the combined tracking and regu- 
lator problem discussed in Section 3.2.3. For this example we used the following state 
space equations for the augmented system: 

z = F(z)z + G(z)u (5.15) 


with state vector: 



X 


flexible manipulator states, see equation (5.8) 

z = 

X* 

= 

desired dynamics states 


e/ 


J edr = integral error terms 


where: 


e 


The dynamics matrices are: 


#i - Ou 
02 — 02 d 



A(x) 

0 

0 


B(x) 

F(z) = 

0 

A d 

0 

G{ z) = 

0 


c 

-Cd 

0 


0 


(5.16) 


(5-17) 


(5.18) 
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where A(x), B(x) are the same as in equation (5.9), while Ad, the desired dynamics 
model, is given by: 


0 10 0 

-u 2 n 0 0 

A d = 

0 0 0 1 

0 0 -u 2 n —2(u>n 

with: 

C = 1-1 


(5.19) 


(5.20) 

(5.21) 


To compute the control input for both the continuous Ricatti design method, and the 
long-term look-ahead controller, we used a cost function of the form: 


J = z T (tj)Hz(t f ) + [ f (ejQjei + e T Q e e + x T Q x x + p 2 u T u)dT (5.22) 

Jt 0 

= f z T (tj)Hz{t f ) + f f (z T Qz + /> 2 u T u)dr (5.23) 

Jt a 

(5.24) 


In both cases we also used the same parameters QuQe^Qx and p. The numerical 
values were: 


Qi 


Qe 

Qx 

P 


5000 0 

0 100 

1 x 10 6 0 

0 1 x 10 5 


diag( 

1 x 10 


0 100 0 10 0 100 0 10 


-4 


(5.25) 

(5.26) 


(5.27) 

(5.28) 
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In order to accentuate the nonlinear nature of the system and demonstrate the capa- 
bility of the control methods used, we chose a reference trajectory which would cause 
the system to perform a maneuver which we do not expect to be executed in reality. 
The desired maneuver (see figure 5-2) starts with the manipulator at a stationary 
position with the arm stretched out (figure 5-2 A). Then the outboard link rotates 
through 360° while the inner link rotates through 90° (figure 5-2 B) until we reach a 
terminal condition where the manipulator is again fully stretched out but now in a 
direction 90° relative to its original orientation (figure 5-2 C). 


Controller Based on Linearized Dynamics 

To demonstrate the nonlinear nature of the problem we first examine the closed 
loop response resulting from a linear design. Using the cost function and controller 
structure described above we found a feedback law by solving an infinite time {tf 
oc) linear quadratic regulator problem using a linearized dynamics model for the 
system. Figure 5-3 shows the resulting unstable behavior. 


Continuous Ricatti Design 

Figure 5-4 shows the well behaved responses we get when we use a continuous Ricatti 
design method to control the flexible manipulator. The control input was computed 
using the system dynamics model and cost function described above. In order to 
reduce the amount of computation we used the method discussed in Section 3.5, i.e. 
we used the feedback gain: 

A'(z) = y 2 R~ 1 G T (z)P old (5.29) 
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Start Condition 




Final Condition 

i 


Figure 5-2: Flexible Manipulator Maneuver 
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Figure 5-3: Two Link Flexible Manipulator Responses - Linearized Design 
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where P 0 i d was found by re-solving the steady state Ricatti equation whenever the 
residual: 


A(z) = Q + PouF{ z) + F T (z)P old - j 2 P old G(z)R- 1 G T (z)P o!d (5.30) 


became large enough that: 

^ > 0.01 (5.31) 

’ M||2 

We see that in this case the continuous Ricatti design controller is able to cope with 
the nonlinear nature of the dynamics. 



Long-term Look-ahead Control 

Finally we applied the long-term look-ahead control method to this example. The 
control input was found as described in Section 4.5, i.e. 

• The expected trajectory was estimated by computing the response of the linear 
time invariant system: 

2 = ^F(zo) - ^G(z 0 )i?- 1 G T (z 0 )Fo^ z (5.32) 

where Zo was the initial condition of the system. 

• The feedback gain was found by integrating the time-varying Ricatti equation: 

_appt I l = g + p{T t )F(i( T )) + FT(z(T))P( T ,t,) (5.33) 

or 

-\p(T,t,)G(i(T))R-'G T (i(T))P(T,tj) (5.34) 

p 

backwards in time. The terminal condition used as a starting point was: 

P{tf, tf) — H — P (5.35) 
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Responses: 2 Link Arm - Continuous Ricatti Design 





Figure 5-4: Two Link Flexible Manipulator Responses - Continuous Ricatti Design 
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where P was found by solving the steady state Ricatti equation: 

0 = Q + PF{ 0) + F t (0)P - \PG(0)R- 1 G T (0)P (5.36) 

P 

Using this value for H ensured a smooth transition to a controller appropriate 
for t > tf, since 

— for t > tf the system will be close to the terminal equilibrium point, 

— and the feedback gain: 

K = ~-R~ 1 G t (0)H (5.37) 

r 

will be appropriate for the system dynamics linearized around the equilib- 
rium point. 


5.4 Summary 

In this Chapter we used both the continuous Ricatti design method of Chapter 4, 
and the long-term look-ahead control method of Section 3 to obtain controllers for a 
complex nonlinear system, i.e. a two-link flexible manipulator. 
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Responses: 2 Link Arm - Long Term Look-ahead 





Figure 5-5: Two Link Flexible Manipulator Responses - Long-term Look-ahead 
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Chapter 6 


Conclusions 


In this thesis we examined methods to synthesize controllers for nonlinear systems. 
Our approach was to exploit the fact that nonlinear systems of the form: 

x = f(x) + G(x)u (6.1) 

could be expressed in quasilinear form, viz: 

x = T(x)x + .B(x) u (6.2) 

provided that f is continuously differentiable [63]. The control methods we developed 
were generic since a wide range of systems could be expressed in the quasilinear form of 
equation 6.2. Furthermore the methods allow the designer to conveniently modify the 
system’s response by adjusting parameters of a cost function. The methods however 
do not guarantee global closed loop stability. 

We developed two classes of control methods: 

• zero-look-ahead control , i.e. control laws that determine the control input based 
only on the current values of T(x), B(x). 
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• controllers with look-ahead, i.e. control laws that use estimates of the future 
behavior of .4(x), B(x) to determine the control input. 

The first type of control law was found by continuously solving a matrix Ricatti 
equation as the system progressed along a state trajectory. We found that by ap- 
propriately adjusting the cost function the method yielded useful control laws. Best 
results were obtained when we cast the problem into a combined tracking/regulation 
problem (see Section 3.2.3) and used very small weights on the control input. With 
regard to stability we found that for this method we could guarantee local stability 
only, but motivated our expectation that the system would be stable for a larger 
range of initial conditions. 

The controllers of the second type used the similarity between quasilinear systems 
and linear time varying systems to find approximate solutions to optimal control type 
problems. Insofar as the control inputs determined in this way were actually optimal, 
the closed loop systems would be stable. 

Our main conclusion is that the methods we developed provide a useful alternative 
which control engineers can use to obtain control laws for systems that exhibit signif- 
icant nonlinearity. However, it also became clear to us that trying to develop control 
methods that are "‘generic'’ , while at the same time guaranteeing global stability pro- 
vides a significant challenge. A more useful approach might be to study nonlinear 
problems on an “individual” basis so as to exploit the “individual physics” of each 
problem to develop control laws (see e.g. [56]). For instance, it may be possible to 
obtain a controller that performs well for the two-link flexible manipulator example 
we considered, by using a feedback law that is “dissipative”. This can be done using 
control torques that are proportional to the joint angles and angular rates. The effect 
would be the equivalent of having a spring and dashpot attached to each joint. Since 
the spring and dashpot combination would be dissipating energy while the system is 
in motion, and would not be adding energy to the system, we would expect the closed 
loop system to be stable and eventually settle. 
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6.1 Future Work 


A key issue regarding any control system is closed loop stability. Ideally a design 
procedure will guarantee (a priori) that the closed loop system will be stable. If this 
is not the case, the stability of the closed loop system has to be determined as part of 
an iterative design cycle. In this context Zubov’s method (see [21] and Section 3.4.2) 
provides a powerful method to assess the size of the stability domain for a given 
design. The results we had show promise, but further work is required to better 
understand the numerical conditioning of Zubov’s equation 3.172. 

In Section 2.3.1 we noted that the system: 

x = A(x)x (6.3) 

would be stable provided that the eigenvalues of A(x) were stable, and that the rate 
of change of the eigenstructure was small enough. This raises the possibility of finding 
stabilizing controllers for systems: 

x = A(x)x + B(x) u (6.4) 

by using an eigenstructure assignment strategy [], which has the aim to keep the eigen- 
structure as constant as possible. Note that the eigenstructure assignment strategy 
also provides the designer with the capability to adjust system responses, see e.g. []. 

Theorem 2.3.1 provides another avenue for obtaining stabilizing controllers for quasi- 
linear systems — the goal would be to obtain a closed loop system matrix with diag- 
onal entries that are negative and dominate the off-diagonal elements. For example, 
we may assume that the control input is given by a feedback law: 

u = -A'(x)x (6.5) 
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This will result in the closed loop system: 


x - A ct (x)x 

= (A(x) - B(x)K(x))x 


( 6 . 6 ) 

(6.7) 


A design procedure could then be: 


• Find A' which minimizes ||£||, where: 


E= A- BK- 


di 0 

0 d n 


( 6 . 8 ) 


This optimization process could be done analytically, resulting in expressions 
for the optimum values of K{d{) and E(di). 


• The next step, which could be done numerically, is to minimize E(d ,) so as 
to make the diagonal elements of the closed loop system matrix as negative as 
possible. Theorem 2.3.1 gives the precise conditions needed for this procedure 
to result in a stable system. 


Finally we note that the optimal control approach has desirable features such being 
“generic” and providing guaranteed stability and robustness. It also allows the de- 
signer to conveniently adjust system responses via a choice of cost function. These 
attractive features should encourage further research in this direction. To make this 
approach more pratical, it will be necessary to find methods that: 

• efficiently solve optimal control problems and provide the answers in terms of a 
state feedback law (for example by solving the HJB partial differential equation). 

• efficiently store state feedback laws found by numerical means. 


163 



Appendix A 


Modeling and Automatic 
Generation of Equations of Motion 
using Maple 

A.l Overview 

This document describes a set of routines useful for deriving equations of motion for 
planar flexible manipulators using the symbolic mathematics software package Maple. 

The following can be modelled: 

• Planar motion of a flexible manipulator (subject to the base being fixed relative 
to inertial space.) 

• A maximum of 9 links (this restriction can easily be removed) 

• Inertial or relative coordinates can be used (see Section A. 2. 2) 

• Each link can be made rigid or flexible, with the following restrictions: 
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- A Bernoulli-Euler beam model assumed. 

- An assumed modes method is used with user supplied mode shape func- 
tions. 

• Fore-shortening effects (stiffening) can be included. (See Section A. 5) 

• With regard to the joints the following can be modelled. 

- Joint masses can be specified for each joint. 

— The effects of joint inertias are not included. 

— Flexibility in a joint is modelled as a simple torsional spring (the joint can 
be made rigid). 

A. 2 Dynamics Model 

In this section we describe the derivation of the equations of motion that were imple- 
mented in the routines. 

Consider a planar manipulator as seen in Figure A-l with axis systems as described 
in Section A. 2. 2 

We can derive the equations of motion using Lagrange’s equations, which are: 

d_ dL _ 31 
dt dqi dq x ~ Ql 


where: 
<h = 

ith generalised coordinate 

Q> = 

ith generalised force 

T = 

kinetic energy 

V = 

potential energy 

L = 

T - V 
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Figure A- 1 : Coordinate 
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Note that if we can find an expression for the kinetic and potential energy in the 
form: 

T = q T Hq 

V = q r Aq 

where: 

q T = [?1 ?2 ■ ■ • ?n] is a vector of generalized coordinates 
H = a generalised inertia matrix 
K = a generalised stiffness matrix 

Lagranges equations can be expressed as: 

Hq + Hq - VT + Kq = Q 

where: 

H = j x (H( q)) is the time derivative of the generalised inertia matrix 
Q T = [Qi Q2 • • Qn] is a vector of generalised forces 

V7T 1 \9T_ dT_ 3T 1 

1^(71 ’ 9^2 ’ * * ‘ ’ d<ln -I 

In the next few sections we will find expressions for the kinetic and potential energy 
of a flexible manipulator. 

A. 2.1 Notation 

In general symbols will have the following meaning: 

• A subscript, or superscript, i, will associate a variable with the ith link. 

• N denotes the number of links of the manipulator. 

• Li denotes the nominal length of link i. 
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A. 2. 2 Axis systems 


Depending on whether we use inertial or relative coordinates the following sets of axis 
systems will be used (see also Figures A-l and A-2). 


Axis Systems for Relative Coordinates 

We have: 

• An inertial axis system Xo,yo fixed at the base of the arm. 

• A rotor axis system x ri , y r , for each link. This axis system is associated with 
the rotor of the torque motor at the root of the ith link (see Figures A-l and 
A-‘2). This axis system is only of importance if we are including flexibility in 
the joints (see also section A. 2. 4 ). The orientation of the rotor axis system can 
be specified by an an angle, 0;, relative to the line tangent to the tip of the 
previous link. 

• A root axis system x t , y , for each link associated with the root of the link. 
In this case the orientation of the axis system is measured relative to the line 
tangent to the tip of the previous link. The angle 0; for this case would typically 
be the angle measured by a joint angle sensor between the previous link and 
the output from a gearbox. The deflection of a link is measured relative to the 
root axis system. 

Relative coordinates would typically be used when: 

• We want state variables which are directly measured (i.e. joint angles). 

• A direct drive arm is being modelled. 

However relative coordinates usually result in more complex equations. 
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Axis Systems for Inertial Coordinates 


We have: 

• An inertial axis system x 0 .y 0 fixed at the base of the arm as before. 

• A rotor axis system x r< , y f| for each link, similar to the case for relative coordi- 
nates. The difference is that now the orientation of the axis system is measured 
relative to an inertial reference (see Figure A-l). 

• A root axis system x, , y« for each link associated with the root of the link. In 
this case the angle 9 t specifies the orientation of the root axis system relative 
to the inertial axis system. The deflection of a link is measured relative to the 
root axis system. 

Inertial coordinates would typically be used when: 

• An indirect drive manipulator is being modelled, i.e. the joints are driven via 
some drive mechanism at the base of the manipulator. 

• We want simpler equations of motion. 

However using inertial coordinates results in equations of motion where the state 
variables are not directly measured. 


Notes 


The deflection of the link is measured relative to the root axis system. By suitable 
choice of the mode shape functions, the x, axis could be, e.g.: 


• tangent to the root of the link, e.g., by using clamped-free cantilever mode 
shapes. 
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• such that the extension of the x t axis joins the root and tip of the link, e.g., by 
using mode shapes of the form: 


x, 


(j)(xi) = sin(nir — ) 

Li 


n = ±1, 2,3... 


where: 

L t = nominal length of the link 
x, = distance measured along x, 

Note that when joint flexibility is taken into account the twist angle between the 
“input and output of the gearbox” will be given by: 

^twist — 

for both inertial and relative coordinates. 


A. 2. 3 Kinetic energy 

Kinetic Energy of the Link (excluding root and tip masses) 

Consider an arbitrary link i and an arbitrary “slice” of this link. The “slice”is located 
at a distance x t along the x, axis when the link is in the undeflected state. We will 
first find an expression for a position vector from the inertial axis system to this 
slice and then take time derivatives etc., to find the kinetic energy associated with 
the slice. The kinetic energy of each link can then be found by integrating over the 
length of the link. 

To keep track of the 'erent axis systems we will use the following notation. A 
vector, R , will be written as: 


170 



RELATIVE CO-ORDINATES 



representing gearbox stiffness 


Figure A-2: Model of joint at the root of link i 


Link (l) 
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(fl) w 

The superscript (i) denotes that the vector is measured relative to axis system x,, y t , 
i.e. when we are dealing with a postion vector, the vector is measured from the origin 
of the axis system indicated by the superscript. We will also, unless stated otherwise, 
assume that the vector is coordinatised in its natural axis system, i.e. if we write the 
vector (R(x t -))^ as a column of numbers, the numbers will give the components of 
(«(*,•))"» in the x,-,y ; frame. 

First find the position vector R(x t )^ from the origin of axis system x,,y ; to 
“slice”x;: When link i deflects, the slice x; moves a distance ix(x,, t) parallel to the 

x, axis and a distance v(x t ,t) parallel to the y, axis (directions for positive movement 
are shown in Figure A-l) (the variable t indicates time dependence). Hence: 

(i?(x t ,0) (,) = ( x t + u(x u t))x % + (v(xi,t))yi 

or in matrix form: 

+ u(xj, t) 
v(xut ) 

Since we are using an assumed modes method we have: 

n, 

v(xi,t) = 53 x)qij(t) 

J=1 

where: are the assumed mode shapes for link i. 

qij are the generalized coordinates for the deflection of link i. 
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Foreshortening effects can be taken into account as shown in Section A. 5. When 
foreshortening effects are taken into account we have: 



When we ignore foreshortening we set: 


u(x,, t) = 0 


To transform the vector to inertial axis systems we will use the so-called 

“homogenous rotation matrices” as described in [3]. To do this we first augment our 
position vector to get: 


Xi -I- u(Xi, t ) 


(*(*,-, 0) (<) = 


v{x,) 


1 


(Note that the superscript has a similar meaning as before). The augmented position 
vector to "slice”x; measured from the origin of axis system x,_x,y;_x and coordinatised 
in axis system x,_x,y,_i is then given by: 




where: 


o-i 


c ;- 1 _„«))<-» 

0 0 1 


and 




L,-i + u(T,_x, t) 
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•_! _ COS [Orel,) -sin(0 re (,) 

C i 

sin(0 re ; t ) cos(9 re i t ) 

and 9 reit is the relative angle between axis systems x x ,y, and x t _ x , y t -i- 
We get the following expressions for 9 re i % : 

• For relative coordinates:(see Figure A-l) 

^reli ^ i "h ^i-1 

where a,_i is the slope (in radians) of link i - 1 at its tip (measured relative to 

• For inertial coordinates: 


9rel, ~ ~ &i - 1 


Assuming small deflections we have: 

dv(xi, t) 

*' = ^t I l - 

We can then transform the augmented vector A(x x , t )<') back to the inertial coordinate 
system by repeated application of the relation above. We finally get: 


X(xi,t)W = A?^...A;- 1 X(i i ,t) (i > 

and extract (#(xi,t))< 0) (position vector from inertial axes to “slice”x t ) as the first 
two components of (A(x,))(°h 

Find an expression for the velocity of the “slice”: Small motions of slice x, 

are related to small changes in the generalised coordinates via the Jacobian matrix. 
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i.e.: 


dR(xi,t )(°* = 7)]rfq 


where: 


J(x, , t ) 


8q 


and q is now a vector of generalized coordinates for the manipulator. In the most 
general case where joint flexibility is included we have: 


01 
0 1 
9n 

<712 

<7l3 


q = 


02 

B 2 

<721 

<?22 


0jV 

<7;vi 

?AT2 


Dividing both sides by dt we get an expression for the velocity of slice x, in inertial 
coordinates viz: 
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R{x„ t) (0) = J(xi,t ) q 


Find the kinetic energy of the ith link: We do this by summing the contribu- 

tions of all the slices in the link, i.e., integrating over the length of the link. (It is 
best to think of x t as a label for the specific slice on the ith link - thus we integrate 
from 0 to Li). The kinetic energy of link i is given by: 

Ti = \ ( pA l q r J T {x i ,t)J(x i ,t)qdx t 
2 Jo 


where: 

pAi is the mass per unit length of link i 
Li is the length of link i 

J(x„t) is the Jacobian of the position vector R(x i: f) (0) indicated above. 


Kinetic Energy of a Point Mass at the tip and root of a link 


Using the same notation as above we see that motion of a point mass at the tip of 
link i can be found by setting: 

x, = L { 


and 


pAdxi = ttil , 


where is the point mass at the tip of link i. 

So. using the same equations as in the previous section, the expression for the kinetic 
energy of a point mass at the tip of link i is: 
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T tipi = ^m i ,q T J T (I 1 )J(I t )q 


Similarly, the kinetic energy of a point mass at the root of link i is: 


Troot, = ^mo,q r J r (I t _ 1 )J(Z:,_ 1 )q 


Total kinetic energy 

Adding the contributions of each link, the total kinetic energy for the system is given 
by: 

AT 

t = £(r t + r tip , + r rooti ) 

1=1 

We have neglected: 

• The effects of rotary inertia of the beam slices 

• The effects of rotary inertia at the joints 


A. 2. 4 Potential Energy 

Potential Energy Associated with the bending of a link 


For a simple Bernoulli-Euler beam model the potential energy associated with the 
bending of link i is given by: 




1 

2 



dx{ 
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where EI t is the bending stiffness of link i. 

Again, using the assumed modes approach we get: 


Vi = \ Eh •)*;(*) j dxi 


which can be rewritten as: 


where: 


Vi = q?\f L ' EIi[V\xi,t)\dxi 
Jo 


~ T __ 
- 


[$"(*,, <)] = 


A 6i 

<7.1 • • • 

9in, J 

On 

VkP* 

... 





C.^2 

C.C, 


Potential Energy Associated with joint flexibillity 

For joint i (joint i appears at the root of link i) we find: 


V jointl = - Oi) 3 

where k x is the joint stiffness of the ith joint. Note that this relation holds for both 
inertial and relative coordinates. 
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Total Potential Energy 


The potential energy for link i and joint i using relative coordinates can then be 
written in matrix form as: 


T Tjojntj — A; <7; 


where: 


K 


k{ — ki 0 

— k{ k l 0 


0 0 


Hr Ei.^dx, 


0 

0 


0 0 ft /„ 1 ’£/,C,Cyx 


Then collect the total Potential energy for the system as: 



AT 

0 . 

..0 

0 . 

..0 . 

.. 0 

...0 

1 T 

Vtot = ^q r 

0 ., 

..0 

AT 

0. 

..0 . 

.. 0 

...0 


0 . . 

,.0 

0 .. 

.0 

0 .. 

o 

o 

.0 

ATv 


Notes: 


• The effects of gravity have been neglected. 

• The effects of shear deformation have been neglected. 


A. 2. 5 Equations of Motion 

Li sing the expressions for kinetic and potent ial energy we find the equations of motion 
are as indicated before: 
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ff(q)q+ ffq-VT + A'q = Q 


It is often convenient to collect the second and third terms in the equation above as: 


C(q,q) = tfq-VT 


A. 3 Using the Routines 

The routines calculate the following terms in Lagrange’s equations for a flexible ma- 
nipulator: 

H( q), ffq, VT, K 

as well as: 

C( q, q) 

The characteristics of each manipulator, for which we want to get equations of motion, 
are specified by filling in the details in a template file such as template, ffl.pl (see 
Section A. 6). Note the following conventions which must be adhered to: 

• Upper and lower case is important. 

• The names of the generalised coordinates must be specified analogous to those 
shown in “template. mpl’\ i.e.: 

- For the rotor angles we must use: 

PSI1, PSI2, 

where the digits specify the link numbers. If joint flexibility is not modelled 
the rotor angles can be omitted totally. 
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- For the angles at the roots of the links we use: 

TH1, TH2, 

— For the deflection coordinates we must use: 

Qll, Q12, .... Q21 , Q22 , 

where the first digit in each name specifies the link number, and the second 
digit in each deflection coordinate specifies the mode number for the given 
link. 

• The names of the mode shapes must be of the form: 

phill, phil2, ,phi21 , phi22, 

where the first digit in each name specifies the link number, and the second 
digit in each shape specifies the mode number for the given link. 

• The names of the mode shapes must be in lower case as shown. 

• Mode shape functions must be specified in Maple operator form, e.g.: 
phill:= < sm(Pi*x/Ll) I x > ; 

To run the routines one must first start Maple with a “clean slate” (e.g. by quitting 
Maple and then starting it again). The user should then specify the file which contains 
the specifications of the system for which the the equations of motion are required. 
This is done by typing, e.g.: 

> modelname := J 'template . mpl ( * ; 

(note that the user can create a specification file with an arbitrary name as long as 
it contains the relevant data in the same form as shown in “template. mpl”). 

To run the routines the user then types: 

> read ‘ armmain . mpl f ; 

This causes the module armmain. mpl to be executed, which reads the data from 
template. mpl and does minor error checking. It also defines some variables, i.e. 
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iner, stiff and lagran which associates the correct file names with the modules which 
calculate the different terms in Lagranges equations. The user can modify arm- 
main. mpl to reflect the proper file names for the user’s site. 

Next the user types: 

> read iner; 

This module calculates an initial form of the inertia matrix H( q). It generally takes 
the longest to execute. 

After the iner module has completed the user types: 

> read stiff; 

This module calculates the stiffness matrix. 

Finally the user types: 

> read lagaran; 

This module calculates the final forms of the terms in Lagrange’s equations. 

Throughout the session the various modules write messages to the screen, indicating 
their progress to the user. Section A. 7 gives an annotated version of a session where 
the routines were used. 


A. 4 Description of routines used 


The main modules are: 
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TEMPLATE. MPL (or equivalent) : Input data specifying the manipulator 

characteristics. The user has to 
supply this. 

armmain. mpl : The main routine for the package 

iner? .mpl : Calculates the inertia matrix 

stiff? .mpl : Calculates the stiffness matrix 

lagran.mpl : Calculates the matrices 

tf(q) 

vr 

C(q, q) 

in Lagranges equations 

Note: 

The “?” denotes digits related to version numbers of the routines. These may change 
with time. 

Subroutine modules are: 

armmsub.mpl : subroutines used by armmain. mpl 

inersub?.mpl : subroutines used by iner? .mpl 

stiffsub.mpl : subroutines used by stiff ? .mpl 


A. 5 Foreshortening 

In this section we use a simple geometric argument from [47] to determine the stiff- 
ening terms in the equations of motion. Alternative derivations are given in [54]. 

Assume that the length of each link, measured along the curved path it makes in 
space, stays fixed, with value, Z,. Thus the projection of the vector from the origin 
of the axes x, y, to the tip of link i will, in general not have an x, axis projection of 

Li. 
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We can find the expression for the motion of the slice labelled x t , in the x, direction, 
by considering an arbitrary element of the link with length As (see figure A-3). For 
the moment we will only consider the x, axis motion of the tip of this element, relative 
to its own root. The motion of the root of the element, will be accumulated by the 
small elements to the “left” of As. We see that the tip of the element moves locally 
inwards (due to the local slope), by a distance Au. Let the projection of the element, 
As, on the x, axis be Aar. Then, by Pythagoras’ theorem and, using the binomial 
expansion for the square root we get: 


As = (Ax 2 + Aw 2 )2 (A.l) 

= (1 + ( A - 2 ) 

- (1 + (A.3) 


so that: 


Au = As — Ax 


(A.4) 


-(-? 
2 W 


(A. 5) 


The cumulative effect of such inward motions is found by integration. This results in 


the expression: 


u(xiJ) - 



ds{ 
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Figure A-3: Foreshortening 

A. 6 Example of “temp late. mpl” 

tt AUTHOR: J.A.Coetsee 

# 12/10/90: Does integrations late as possible 

# Maple v4.3.1 (IBM-PC) 

« 

# 

g**************** ****************************************** 

# FILL OUT THE FOLLOWING DATA TEMPLATE FOR EACH LINK * 

g******** ********************* ***************************** 


# ! ! Use the I PRESCRIBED CAPITAL I letters for generalised co-ordinates 
« I LOWERCASE I " " mode shapes 


185 




# 


NOTE: 


# ASSUMES mode shapes are orthogonal per link 

# (see procedure "orthorel" in file "linksub? .mpl") 

# 





# 



# DESCRIPTION: 

2 link flexible arm 


« 

Inertial Co-ordinates 


# 

Sinusoid mode shapes 


# 



# DATE : 

11/16/90 


# 






COORDS:* inertial; 

# Co-ordinate system used: 

# "relative":- 

link orientation THi 


* 

measured relative to 


# 

tip of previous link 


# "inertial":- 

link orientation THi 


# 

measured relative to 


# 

inertial frame 


FORESHORTEN : =no 

1 

# foreshortening effect 

# 



NL : =2 ; 

# 


# Number of links 


# Data template for LINK 1: 
nl : = 1; 

qstarl := [TH1 ,Q11] ; 

LI := L; 
mO 1 : = 0 ; 
mLl : = mO ; 


# Number of modes in the link 

# List of generalised co-ords. link 1 

# Link length 

# Root mass 

# Tip mass 
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kl:= infinity; 

Ell : =EI ; 
rhoAl:= rhoA; 

phill:= <-sin(Pi*x/L) I x >; 
phil2:= <-sin(2*Pi*x/L) I x >; 

# 

# Data template for LINK 2: 
n2 : = 1 ; 

qstar2 : = [TH2 , Q21] ; 

L2 := 12; 
m02:= 0; 
mL2:= m; 
k2:= infinity; 

EI2 : = El; 
rhoA2:= rhoA; 

phi21:= <-sin(Pi*x/12) | x >; 
phi22:= <-sin(2*Pi*x/12) | x >; 

*************************************** 


# Joint stiffness 

# Beam stiffness 

# Beam density per unit length 

# Mode shape 1 

# Mode shape 2 


# Number of modes in the link 

# List of generalised co-ords. link 2 

# Link length 

# Root mass 

# Tip mass 

# Joint stiffness 

# Beam stiffness 

# Beam density per unit length 

# Mode shape 1 

# Mode shape 2 

*********************** 


A. 7 Annotated Example of Session 

IWI 

- _ I \ I I / I _ • Copyright (c) 1989 by the University of Waterloo 
\ MAPLE / Version 4 . 3 . 1VM/80387 - February 1990 
< > Fl-Help F2-Search F3-Quit F4-Shell 
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# First give the name of the file specifying the arm for which we want 

# equations of motion 

>modelname:= 'template. mpl ' ; 

modelname := template.mpl 

# Now run the main routine which sets everything up 

# (The routines inform the user of their progress ) 

>read ‘ armmain .mpl ' ; 


Main program - Does initialization 


Warning: new definition for trace 

printlevel := 1 

q := □ 

q := [TH1 , Qll] 

q := [TH1, Qll, TH2 , Q21] 

\ 

i : = i 

q := array (1 .. 4, [TH1 ,qil ,TH2,Q2l] ) 
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qstarsize := 2 
qstarsize := 2 
inertia := iner6.mpl 
stiffness := stiff. mpl 
lagrange := lagran.mpl 

® Then run the file which calculates the initial form of the generalised 
# inertia matrix 

>read inertia; 

I Calculating the inertia matrix | 


Busy with link , 1 


doing v(x) 

bytes used=452504, alloc=204800 , time=3.430 
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doing u(x) 


doing A 


doing RO 


doing Jacobian 


doing JJ 


doing initial integration of JJ 


time for integration of Jacobian is:, .170 


busy with Hiroot 


time for simplifying Hiroot is, .170 


busy with Hitip 


adding up the components of the H-matrix 


bytes used=854452 , alloc = 360448 , time = 6.840 


Busy with link , 2 


doing v(x) 
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doing u(x) 


doing A 


doing RO 


doing Jacobian 


doing JJ 


bytes used=1254688, alloc=483328, time=10.190 

doing initial integration of JJ 


time for integration of Jacobian is:, .110 


busy with Hiroot 


time for simplifying Hiroot is, 0 


busy with Hitip 


adding up the components of the H-matrix 


Starting final simplification of H matrix 


- getting constansts outside integral signs 


bytes used=1655236, alloc=565248, time=13.650 

- applying orthogonality conditions 
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- doing final integrations with the actual mode-shapes 


bytes used=2058144, alloc=565248, time=18.760 
bytes used=2458276 , alloc=565248, time=24.470 
bytes used=2858508, alloc=565248 , time=30.180 
bytes used=3259148, alloc = 565248, time = 34.520 
bytes used=3659364, alloc=565248, time=38.480 

-doing final collection of terms 

Time for last few simplifications of H is , 24.280 

printlevel := 1 


# The first routine has just completed 

# At this point the user can examine values etc. e.g. 


>H[l,ll ; 


2 3 2 2 2 

1/2 rhoA Qll L + 1/3 rhoA L + mO L + rhoA L 12 + m L 


# Now run the file which calculates the stiffness matrix 


>read stiffness; 
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1 Calculating the 

stiffness matrix | 


Busy with 

link , 1 


bytes used=4059568 , alloc=606208 , time=41.860 


[ o 

0 0 

0 ] 

[ 


] 

[ 

4 

] 

[ 

Pi 

] 

[ o 

1/2 — 0 

0 ] 

K-matrix , [ 

3 

] 

[ 

L 

] 

[ 


] 

[ o 

0 0 

0 ] 

[ 


] 

[ o 

0 0 

0 ] 

Busy with 

link , 2 


[ 0 

0 0 

0 ] 
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6-3 




[ ] 

C 4 ] 

[ Pi ] 

[ 0 1/2 — 0 0 ] 
C 3 ] 

[ L ] 

K-matrix , [ 3 

[ 3 

C 4 ] 

[ Pi 3 

[ 0 0 0 1/2 — ] 

[ 3 ] 

[ 12 ] 


printlevel := 1 

# The second routine has just completed 

# Then calculate the terms in Lagrange's equations 
>read lagrange; 


Calculating the terms in Lagrange’s equations 


bytes used=4462016 , alloc=606208 , time=45.760 


Simplifying the "Coriolis" Terms 
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bytes used=4862464, alloc=614400 , time=50.370 
bytes used=5262596, alloc=614400, time=54.930 
bytes used=5662784, alloc=614400 , time=59.760 

Collecting terms in a useful order 

bytes used=6063704, alloc=638976, time=64.210 

Combining terms like: sin(a) *cos (b) + sin(b)*cos(a) > 

bytes used=6463844, alloc=679936 , time=67.890 

printlevel := 1 


# The terms now show time dependence e.g. 

>Coriolis [2] ; 

/ d \2 

- 1/2 L I thl(t)| rhoA qll(t) 

\ dt / 

>H[l , 1] ; 

2 3 2 2 

1/2 rhoA qll(t) L + 1/3 rhoA L + mO L + rhoA L 12 

# The results can be saved to a file if need be e.g. 

>save HH, Coriolis, K, modelname, 'results .out ‘ ; 

# To end the session 
> quit ; 


sin(a+b) 


2 

+ m L 
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Appendix B 


Quasilinear Equations of motion 
for Double Pendulum Example 


The equations of motion for the double pendulum of Example 2.2.1 are: 


tf(q)q + C(q,q)q+ Kq + G grav {q)q = 0 (B.l) 

where: 


H{ 1,1) 

= mi l\ + m 2 l\ + 2 m 2 h h sin(0i) sin(#i + 0 2 ) + m 2 l 2 

(B.2) 


+2m 2 l 2 l\ cos(^i) cos(^i + 0 2 ) 

(B.3) 

H{ 1,2) 

= m 2 1 2 1\ sin(^i) sin(^i + 0 2 ) + m 2 1\ 

(B.4) 


+m 2 l 2 li cos(^i) cos(^j + 6 2 ) 

(B.5) 

H{ 2,1) 

= H{ 1,2) 

(B.6) 

H( 2,2) 

= m 2 l'l 

(B.7) 

C( U) 

— 2 m 2 1 2 U 9 2 sin(^i) cos(#i + 0 2 ) 

(B.8) 


—2m 2 l 2 li9 2 cos(6\) sin(^i + 0 2 ) 

(B-9) 

C( 1,2) 

= m 2 l 2 1\ d 2 sin^) cos(0! + 6 2 ) 

(B.10) 



PMpGSCW*** PAGE BLANK NOT FH.MtD 
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— m 2 1-2 h &2 cos (0i ) sin(#i + 0 2 ) 
(7(2,1) = — 2 0im 2 / 2 /i sin(0i) cos(0i + 02 ) 

+2 0i m 2 / 2 / 1 cos(0i ) sin(0i + 0 2 ) 
C(2,2) = — 0i m 2 / 2 /i sin(0i) cos(0i + 0 2 ) + 
0i m 2 / 2 /i cos(0i) sin(0i + 0 2 ) 


A' = 


A'i 0 

0 K 2 


G 


grav 


V 1 IS. 

9i 


(9i+*2) 


m 2 gl 2 sinf^i +6 2 ) 

(Oi+e 2 ) 


m 2 gl 2 sin(^i +O2 ) 
(#i +# 2 ) 

m 2 gl 2 sin(fli+fl 2 ) 
(#1 +# 2 ) 


(B.ll) 

(B.12) 

(B.13) 

(B.14) 

(B.15) 


( B. 16) 
(B-17) 


198 



Bibliography 


[1] W.F. Ames. Numerical Methods for Partial Differential Equations. Academic 
Press, 1992. 

[2] B.D.O. Anderson. Stability results for optimal systems. Electron. Lett., 5:545, 
Oct 1969. 

[3] H. Asada and J. Slotine. Robot Analysis and Control. Wiley, 1986. 

[4] M. Athans. Class notes for multivariable control systems. 

[5] R.H. Bartels and G.W. Stewart. Algorithm 432: Solution of the matrix equation 
AX + X B = C . Commun. ACM , 15:820-826, 1972. 

[6] E. Bayo, P. Papadopoulos, J. Stubbe, and M.A. Serna. Inverse dynamics and 
kinematics of multi-link elastic robots: An iterative frequency domain approach. 
Int. J. of Robotics Research , 8(6):49-62, Dec 1989. 

[7] R.W. Brockett. Feedback invariants for nonlinear systems. IFAC Congress , 
(6):1 1 15—1 120, 1978. 

[8] A.E. Bryson and Y.C. Ho. Applied Optimal Control. Halstead Press, 1975. 

[9] J.R. Bunch and D.J. Rose, (eds.). Sparse Matrix Computations. Academic Press, 
1976. 

[10] R.H. Cannon, Jr. and Schmitz E. Initial experiments on the end-point control 
of a flexible one-link robot. Int. J.of Robotics Research, 3(3):62-75, 1984. 


199 



[11] A.R. Collar. Some notes on Jahn s method for the improvement of approximate 
latent roots and vectors of a square matrix. Quart. J. Mech. Appl. Math , 1:145— 
148, 1948. 

[12] W.A. Coppel. Stability and Asymptotic Behavior of Differential Equations. 
Heath, 1965. 

[13] G. Dahlquist. Stability and error bounds in the numerical integration of ordinary 
differential equations. Trans. Roy. Inst. Tech. (Sweden), 130, 1959. 

[14] C.A. Desoer and H. Haneda. The measure of a matrix as a tool to analyze com- 
puter algorithms for circuit analysis. I.E.E.E. Transactions on Circuit Theory. 
19(5):480— 486, 1972. 

[15] P. Dyer and S.R. McReynolds. The Computation and Theory of Optimal Control. 
Academic Press, 1970. 

[16] B.A. Frances and G. Zames. On H optimal sensitivity theory for SISO feedback 
systems. I.E.E.E. Transactions on Automatic Control , 29:9-16, 1984. 

[17] J.H. Ginsberg. Advanced Engineering Dynamics. Harper and Row, 1988. 

[18] G.H. Golub and C.F. Van Loan. Matrix Computations. Johns Hopkins University 
Press, 1983. 

[19] G.C. Goodwin and K.S. Sin. Adaptive Filtering Prediction and Control. Prentice- 
Hall, 1984. 

[20] J.D. Graham, R. Ravindran, and K. Knapp. Space manipulators - present capa- 
bility and future potential. In AIAA/NASA Conference on Advanced Technology 
for Future Space Systems , pages 243-253, 1979. 

[21] W. Hahn. Stability of Motion. Springer- Verlag, 1963. 

[22] R.M. Hirschorn. Invertibility of multivariable nonlinear control systems. IEEE 
Trans, on Aut. Control , AC-24:855-865, 1979. 


200 



[23] L.R. Hunt, R. Su, and G. Meyer. Design for multi-input nonlinear systems. In 
Differential Geometric Control Theory. Birkhauser, 1983. 

[24] A. Isidori. Nonlinear Control Systems. Springer- Verlag, 1989. 

[25] B. Jakubczvk and W. Respondek. On linearization of control systems. Bul- 
letin de L ’Academie Polonaise des Sciences, Serie des sciences mathematiques , 
XXVIII:517-522, 1980. 

[26] T. Kailath. Linear Systems. Prentice-Hall, 1980. 

[27] S. Kenny. Private communications. 

[28] D.E. Kirk. Optimal Control Theory — An Introduction. Prentice-Hall. 1970. 

[29] D.L. Kleinman. On an iterative technique for Ricatti equation computations. 
IEEE Trans, on Automatic Control , AC-13:1 14-115, February 1968. 

[30] P.V. Kokotovic, H.K. Khalil, and J. O’Reilly. Singular Perturbation Methods in 
Control: Analysis and Design. Academic Press, 1986. 

[31] A.J. Krener. On the equivalence of control systems and the linearization of 
nonlinear systems. SIAM J. Control , 11(4):670— 676, 1973. 

[32] H. Kwakernaak and R. Sivan. Linear Optimal Control Systems. Wiley- 
Interscience, 1972. 

[33] H. Kwakernaak and R. Sivan. The maximally achievable accuracy of linear 
optimal regulators and linear optimal filters. IEEE Transactions on Automatic 
Control , AC-17:79-86, 1972. 

[34] W.H. Kwon, A.M. Bruckstein, and D.G. Byun. Receding horizon tracking control 
as a predictive control and its stability properties. Int. J. Control , 50:1807-1824, 
1989. 

[35] W.H. Kwon, A.M. Bruckstein, and T. Kailath. Stabilizing state-feedback design 
via the moving horizon method. Int. J. Control , 37:631-643, 1983. 


201 



[36] W.H. Kwon and A.E. Pearson. A modified quadratic cost problem and feedback 
stabilization of a linear system. IEEE Transactions on Automatic Control , AC- 
22:838-842, 1977. 

[3<] J. La Salle and S. Lefschetz. Stability by Liapunov’s Direct Method. Academic 
Press, 1961. 

[38] A.J. Laub. A Schur method for solving algebraic Ricatti equations. IEEE Trans, 
on Automatic Control , AC-24(6):913-921, December 1979. 

[39] A.G.J. MacFarlane. An eigenvector solution of the optimal linear regulator prob- 
lem. J. Electron. Contr ., 14:643-654, 1963. 

[40] K. Martensson. On the matrix Ricatti equation. Information Sciences , 3:17-49. 
1971. 

[41] D.Q. Mayne and H. Michalska. Receding horizon control of nonlinear systems. 
IEEE Trans, on Automatic Control , AC-35(7):814-824, July 1990. 

[42] 0. Mayr. The Origins of Feedback Control. M.I.T. Press, 1970. 

[43] R.H. Middleton and G.C. Goodwin. Adaptive control of time-varying linear 
systems. I.E.E.E. Transactions on Automatic Control , 33(2):150-155, 1988. 

[44] D. Mitra and C. So Hing. Linear inequalities and p matrices with applications to 
stability of nonlinear systems. In Proc. of the 5th Asilomar Conf. Circuits and 
Systems , pages 179-188, 1971. 

[45] P.J. Moylan and B.D.O. Anderson. Nonlinear regulator theory and an inverse 
optimal control problem. IEEE Trans, on Automatic Control, AC-18(5):460-465, 
October 1973. 

[46] J.R. Munkres. Analysis on Manifolds. Addison- Wesley, 1991. 

[47] Carlos E. Padilla. Nonlinear strain-displacement relations in the dynamics of a 
two-link flexible manipulator. Master s thesis, Department of Aeronautics and 
Astronautics, M.I.T., 1989. 


202 


[48] Carlos E. Padilla. Performance Limits and Robustness Issues in the Control of 
Flexible Link Manipulators. PhD thesis, Department of Aeronautics and Astro- 
nautics, M.I.T., 1992. 

[49] J.-H. Park and H. Asada. Design and analysis of flexible arms for minimum-phase 
endpoint control. American Control Conference , 1990. 

[50] W.H. Press, B.P. Flannery, S.A. Teukolsky, and W.T. Vetterling. Numerical 
Recipes in C. Cambridge University Press, 1988. 

[51] H.H. Rosenbrock. A Lyapunov function with applications to some nonlinear 
physical systems. Automatica. 1:31-53, 1962. 

[52] H.H. Rosenbrock. A Lyapunov function for some naturally ocurring linear ho- 
mogeneous time-dependent equations. Automatica , 1:97-109, 1963. 

[53] J.S. Shamma. Analysis and Design of Gain Scheduled Systems. PhD thesis, 
Massachusetts Institute of Technology, 1988. 

[54] J.C. Simo and L. Vu-Quoq. The role of non-linear theories in transient dynamic 
analysis of flexible structures. Journal of Sound and Vibration , 119(3):487-508, 
1987. 

[55] J-J.E. Slotine. Sliding controller design for nonlinear systems. Int. J. Control , 
40(2), 1984. 

[56] J-J.E. Slotine. Putting the physics into control. IEEE Control Systems Magazine , 
8(6), 1988. 

[57] J-J.E. Slotine and J.A. Coetsee. Adaptive sliding controller synthesis for nonlin- 
ear systems. International Journal of Control , 43(6) : 163 1—1651 , 1986. 

[58] J-J.E. Slotine and W. Li. On the adaptive control of robot manipulators. Int. 
J. Robotics Research , 6(3), 1987. 

[59] J.J. Slotine and W. Li. Applied Nonlinear Control. Prentice-Hall, 1991. 


203 



[60] G. Stein and M. Athans. The LQG/LTR procedure for multivariable feedback 
control design. IEEE Trans, on Automatic Control, AC-32(2):105-1 14, 1987. 

[61] R.F. Stengel. Stochastic Optimal Control. Wiley Interscience, 1986. 

[62] P. Van Deventer. Pseudo-linear control methodology for nonlinear systems. Tech- 
nical report, Massachusetts Institute of Technology, 1993. 

[63] M. Vidyasagar. Nonlinear Systems Analysis 2nd ed. Prentice- Hall, 1993. 

[64] V. I. Zubov: Methods of A. M. Lyapunov and Their Application. NoordhofF, 1964. 


204 



